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1. INTRODUCTION , 


A numerical model for simulation of the global circulation of the 
stratosphere and mesosphere is currently under development at the University 
of Washington. The complete model is a semi-spectral model in which the 
longitudinal dependence is represented by expansion in zonal harmonics 
while the latitude and height dependencies are represented by a finite 
difference grid. Since many of the dynamical processes which occur in the 
stratosphere and mesosphere are the result of interactions between the 
zonal mean flow and planetary waves , it is useful to formulate a model in 
which the zonal mean and wave portions are explicitly separated, as is 
done here. 

The model is based on the primitive equations in the log pressure 
coordinate system as given by Holton (1975). In order to avoid the problems, 
inherent in simulating tropospheric meteorological processes, the lower 
boundary of the model domain is set at the 100 mb level (i.e., near the 
tropopause) and the effects of tropospheric forcing are included in the 
lower boundary condition. The upper boundary is at approximately 96 km, 
and the latitudinal extent is either global or hemispheric. 

In this report we first outline the basic differential equations 
and boundary conditions. We next describe the finite difference equations. 
We then discuss the initial conditions and present a sample ’calculation. 
Finally, the Fortran code is given in the appendix. 



2. BASIC EQUATIONS 


In setting down the basic equations we will make use of the following 
symbols ; 

X longitude 

6 latitude 

z a measure of "height" [= -H Jin (p/p^)] 

H scale height [e RT^/g] 

R gas constant for dry air 

Tg a constant stratospheric mean temperature 
g gravitational acceleration 

p pressure 

p^ a constant reference pressure 
u eastward velocity 

V northward velocity 

w a measure of "vertical velocity" [e dz/dt] 

a basic state temperature [e T^(z)] 
a basic state geopotential [= 'l'^(z)] 

T departure of local temperature from 1 ^( 2 ) 
departure of local geopotential from ^^( 2 ) 

Q angular velocity of earth 

a radius of earth 

J diabatic heating rate per unit mass 

Cp specific heat at constant pressure 

K ratio of gas constant to specific heat at 
constant pressure [e R/c^] 

dx eastward distance increment [e a cos 6 dX] 

dy northward distance increment [= ad6]. 
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The horizontal momenttun equations can then be written in flux form as 


3u 

3t 3x cos^ e 9y 


^ (uv cos^ 0 ) 


H (p wu) - 2Q^7 sxn 0 

p 9z o 
o 




( 2 . 1 ) 


3v . 3 , s . 1 3,2 

•^ + tt* (iiv) + ^ - 5 — (v^ cos 6 ) 

3t 3x cos 6 3y 





tan 0 


+ 2 J^u sin 0 = “ 


3y 


+ D2(v) 


( 2 . 2 ) 


Here, p^ = p^ exp(-z/H), where is the mean density at z = 0 . 

D 2 (v) represent subgrid scale momentum diffusion. Explicit forms for these 
terms will be given in Section 4. 

Using the above notation the hydrostatic approximation and continuity 
equation become 


df RT 
o o 

dz ~ H ’ 


3z H ’ 


(2.3) 


and 



— (v cos 0 ) + — (p w) = 0 
cos 0 3y p 3z o 

o 


(2.4) 


The variables T^ and define a hydrostatically balanced basic state 
which is specified to be the U.S. standard atmosphere. Using (2.3) we can 
write the thermodynamic energy equation for the departure from the basic 
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state as. follows 


1 


3^ 
z 

3t 


+ 


_3_ 

3x 


(u*^) 


cos 0 3y 


(v$ cos 
z 


6 ) 


+ — (p $ w) + wN^ = kJ/H + 1)^(9 ) (2.5) 

p 3z ^o z 2 z 

o 


where 





/dT 
o 

dz 


+ 


kT 
o 

H 

J 


is the buoyancy frequency squared, and we have let denote the subgrid 

scale diffusion. 

The basic state temperature profile is assumed to be in radiative 
equilibrium (see Section 9), so that the horizontal average of the diabatic 
heating will vanish provided that the horizontally averaged temperature 
equals the basic state temperature T^(z). Because of the nonlinearity of 
(2.5) the horizontally averaged temperature need not remain equal to T^(z) 
as the flow evolves in time. However, in practice we find that departures 
of the horizontally averaged total temperature from (z) are at most a few 
degrees so that for practical purposes the horizontally averaged diabatic 
heating remains very small, and J can be regarded as the differential 
heating. 


^Following Holton (1975) we here neglect the small term wkT/H compared to 
wkT^/H. This approximation is necessary if we wish to define available 
potential energy in terms of the temperature variance. 



3. ZONAL HAEMONIC EXPANSION 


The basic equations of the model were given as (2.1), (2.2), (2.4), 
and (2.5). In order to develop the semi-spectral model we expand the 
basic equations in zonal harmonic series by letting 

n =+<» 

f(X,y,z,t) = ^ F^(y,z,t)e^^^ (3.1) 

n = -<» 

where f(X,y,z,t) stands for any field variable and F is the Fourier transform 

n 

of f defined by 


F 

n 


,-z/2H 

2tt 


+ TT 


f(A,y,z,t)e dX 


-IT 


(3.2) 


so that F_^ = F*, where the asterisk denotes the complex conjugate. 

To transform the nonlinear terms in the basic equations we need the 
convolution theorem: 


_L 

2tt 


+ ir 

[f(X)g(X)]e“^"'^ 

-r 


dX 


+ CO 

G F 
m n-m 

m = -«> 



from which we find 


Jl. 

2tt 


F 

<ro 

J 

3X ® 


-inX 

e 


dX 


z/H 


-f- 03 

z 


im G 


n-m m 


- IT 


m = _oo 
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To Fourier transform (2.1), (2.2), (2.4), and (2.5) we define the 
following transform pairs: 


f(A): u V w f kJ/H 

F : U V W ¥ Q 

n n n n n n 


The transformed equations are as follows: 


3U 


n 


3t 


fV = ^ '1' + D„ (U ) 

n a cos 0 n x. n 


+ 00 


- e 


z/2H 


y U U 

L-, [a cos 0 m n- 


■m 


m = 


+ — -I- (U V cos^ 0) + (U W ) 
cos*^ 0 dy m n-m dz m n-m 


(3.3) 


where £. = 1 for n = 0, f. = 2 for n 0. 


9V 


3 ¥ 


■^ + fU = - ^ + D, (V ) 
dt n dy 2 n 


+ CO 
m =— CO 



(U V 
m n-m 


+ V U ) 
m n-m 


+ — ^ ^ (V V cos 0) 
cos 0 dy m n-m 


+ (V ¥ ) + (U U ) 

dz m n-m^ a m n-m 


(3.4) 
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_3_ 

at 



X. 

2H 


n 


+ N^W = 
n 


Q_ + D, 


n 




CO 

^z/ZH ^ 


m = 


r im 

TT 



|;a cos 0 

U 

m 

3z 

2hJ 


+ 



X 

2H^ 


'J' U 
m n-m 


+ 


1 _ 3 _ 

cos 0 3y 


|cos 


6 V 

n-m 







(3.5) 


inU 

n_ 

a cos 0 


— ^ ^ (cos 0 V ) + 
cos 6 dy n Idz 


2H 

J 


W = 0 
n 


(3.6) 


We now severely truncate the wave spectrum by assuming that the flow 
consists of a single wave of wavenumber n = s, and the zonal mean n = 0. 

To exclude all other wave modes we must replace the summations in (3.3) - (3.5) 
by a summation over the two values m = 0 and m = s. 


3.1 The zonal mean equations 

If we set n = 0 in (3.3) - (3.6) and replace ( )^ by ( ) for all field 
variables we obtain the zonal mean equations: 


f - fv = - 

at 


cos 


^ (H V cos^ 0) + ^ cu W) 


+ + D^(U) 


(3.7) 
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(3.8) 


^ ^ (V cos e) + 


cos 6 3y 


3z 


2H 


W = 0 


(3.9) 


a 

5t 



1 + + F ■ 

f ^ ^ 

V cos 0 




9z 

HJ T 

(cos 0 3y 

Sz 

^ 2HJ 


3z 


W 


3z 2H 


+ D, 


^ + -1- + Q 
3z 2H ^ 


( 3 . 10 ) 


where f = 2(2 sin 6 is the Coriolis parameter. Here denotes the convergence 
of the momentiiin flux due to zonally asymmetric motions (e.g., planetary 
waves) while denotes the convergence of the eddy heat flux. 

We have neglected the advection by the mean meridional circulation and 
the eddy momentum flux terms in (3.8) since the mean zonal wind is nearly 
in gradient wind balance. The terms 3V/3t and (V) are also very small 
but must be retained for our method of numerical solution. 

With the aid of (3.9) we can define a mean meridional streamfunction, 

X, by letting 


T7 fi 3X 

W cos U = 

3y 


V cos 0 


3z 2H 


X 


(3.11) 


The X field proves useful in specifying boundary conditions and solving the 
zonal mean component equations. 

The eddy flux convergence terms in (3.7) and (3.10) have the following forms 


^ 


z/2H 


{-a 

(cos'' 


[(U V* + U*V ) cos^ 6] 


2 0 3y '^"s s ' s s 


+ ^ (U W* + U*W ) 
9z s s s s 


(3.12a) 
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and 


= - e 


z/2H 


/_j^Ar 

l^cos 0 3y L 


cos 9 


V + ^! + V* 

s dz 2H s s 


3z ^ 2H 




+ 




9z 2H 


Y + W 
s s 


+ — lY^ 
9z 2Hrs 


:)''s]} 


(3.. 12b) 


3.2 The eddy equations 

When we set n = s in (3.3) - (3.6) and again designate zonal means by 
an overbar rather than the m = 0 subscript, we obtain the eddy equations, 


3U 

_J. _ f V = Y 

0t s a cos 6 s 


z/2H 

e 


■ isU U V p, _ 

^ ^ (U cos 0) 

a cos 6 cos 0 3y 


+ W 

s 



U> + 


”2^3) 


C3.13) 


av 

ir + f's- 


8Y flV U tan 0 isV U 

s _ ^^/2H / s ^ s , + B (V > 

9y la a cos 0] 2s 


(3.14) 


A rfA+ 

9t L(9 z 2hJ^s_ 


+ N^W = Q - e‘ 
s s 


:/2H 


f isU 

'9 1' 

(a cos 0 

9z 2H 




2 _ 3 

Y + V ^ 
s 9y 


9z 2HJ 


Y 


+ B. 


9z 2 h| '^s 


(3.15) 


isU ^ j. f 

(cos 0 V ) + 

a cos 6 cos 0 9y s 


_ 9 _ 

9z ” 2H 


W = 0 
s 


(3.16) 


Here all terms involving advection by the mean meridional circulation, V, W 
have been neglected. 



4. BOUNDARY CONDITIONS 


Conditions for the zonal mean: 

The model can be integrated on either a hemispheric or global domain. 
For global integrations the horizontal boundary conditions are as follows: 

X = U = V = dY/dy = 0 at 0 = ± tt/2 (4.1) 


For hemispheric integrations the same boundary conditions are used at 
0 = 0, TT/2 except that a value of U not equal to zero may be specified 
at 0 = 0. 

Vertical boundary conditions are specified as follows : 

U = U^(y,t) at z = 0 (4.2a) 

where z - 0 designates the lower boundary (i.e., the tropopause level) and 
U is an externally specified mean zonal wind. The boundary mean zonal . 
flow is assumed to be in gradient wind balance. Thus from (2.10) we see 
that at z = 0 

V = 0 (4.2b) 


and 


f . fir + P tan_e ^./2H 
oy a 


(4.2c) 


Using the conditions (4.2a) we can integrate (4.2c) to obtain 'J'(y,t) at 
z = 0, provided that we let the horizontal average of W(y,Q,t) vanish. 

At the upper boundary (z = z^) we assume that the vertical shear of 
the mean zonal wind, the mean meridional wind, and the mean geopotential 
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all vanish . Thus , 



(4. 3a) 



(4.3b) 


and 



Y = 0 


(4.3c) 


Condition (4.3c) of course implies that the zonal mean temperature must 
equal the basic state Tq(z^) at z = z^. 

In addition to these conditions it is clear from (3.7) and (3.10)' that 
boundary conditions are also required for the vertical momentum and heat 
fluxes associated with the mean meridional circulation. We wish to avoid 
specifying W or the fluxes themselves at z = 0. Instead we assume that 
the flixx divergences vanish at the lower boundary: 

^(OW) + f 

However, for simplicity we assume that the fluxes themselves vanish at the 
upper boundary. If in addition we let Q = = 0 at the upper boundary, 

then from (3.10) we have 

W = 0 at z = z^ 


= 0 


(4.4) 


(4.5) 
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Boundary conditions for the eddy equations ; 

The boundary conditions for the eddy motions are analogous to the 
conditions for the zonal mean. However, the case s = 1 must be treated 
separately because does not vanish at the poles for s = 1. Thus, for 
integrations on the global domain we have at 0 = ± tt/2: 

'H = 0 
s 

U = V = 0fors>l (4.6) 

s s 

3u^/ae = av^/ae = o for s = i 

For a hemispheric domain the conditions at 9 = 0 depend on the symmetry 
conditions assumed. If geopotential is symmetric we have 

3'i' /36 = 9U /ae = V = 0 at 0 = 0 (4.7) 

s s s 

If geopotential is antisymmetric we have 

W = U = 9V /30 = 0 at 0 = 0 (4.8) 

s s s 

Conditions at the horizontal boundaries are specified as follows ; 

At the lower boundary a geopotential height perturbation is specified so 
that 


T (y,t) = gh (y,t) at z = 0 (4.9) 

o S 

while at the upper boundary the wave perturbations are assumed to vanish 

T = 0 at z = z„ 
s T 


C4.10) 
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The latter condition requires that we impose strong damping in the layers 
near to prevent spurious reflection of wave energy from the upper 
boundary. Finally, in order to compute and at the upper and lower 
boundaries we assume that the vertical momentum and heat flux divergences 
vanish at the boundaries . 



5. ENERGETICS 


It can be shown that the eddy equations (3.13) - (3.16) are energetically 
consistent with the mean flow equations (3.7) - (3.10). In fact the system 
is governed by a Lorenz type energy cycle which (neglecting the diffusion 
terms) may be written as follows: 


II = (k^ -K>+ (a-K>+ B(K) 


(5.1) 


dt 


= - <A*k)-(a*a) + G + B(A) 


(5.2) 


dK 
s 

dt 


^ - (K • k)+(a • K ) + B(K ) 


(5.3) 


dA 


dt 


® = <A « A^) -<A • K ) + G 


s S' 


(5.4) 


where 



o o 


cos 0 d6 dz 


“Tr/2 

_ ■ f 1 f9V . n ^ 


o o 


<K3-k>e- j 


7T/2 

' We"'® SI 


o o 


+ cos 6 7 ^ (U W* + U*W ) }• d0 dz 
9z s s s s ' 


>} 
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<A.K>= I j 


Tt/2 


Trf 8 . 1 


W 


8z 2H 


T cos 0 d0 dz 


o o 


7t/2 


B(K) = 


{[WC'F + u 2/2) + U(U W* + U*W )] cos e d0} 

s s s s z = 0 


(a • A > 
' s' 


I o 


N2 


fB'l' . Y] / 9 / - L 


[9z 2 hJ'^s 


o o 


+ V* 
s 


9z 




]} 


+ cos 0 “ W* ^ 

9z s 9z 2H 


Y 


+ W 




tt/2 


G E 


E + 

■ 


_s_ 

'3Y Y ' 


4 


N2 

9z 2Hj 


cos 0 d0 dz 


o o 


B(A) 


tt/2 


f j 

’ 1 

r|Y^xl 

^ — 1 
w + — 

'3Y Y' 

f 

w* 

f- - + -1 

J 1 

2N2 

[3z 2 hJ 

N2 

3z 2H 

W y/ 

s 

s 

3z 2H 


Y 


+ W 



^ + 9dei^.„ 


K H 
s 


“ tt/2 


o o 
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1-"/^^ 1 r®'*' 

= J J p- cos 6 dB dE 


o o 


¥/2 


iK ■ kJ = 



• 

w 

3 1 1 

^ + W* 

[- + -1 

’i' 

J 


s 

3z 2H s s 

k / 

9z 2H 

s 


COS 6 d0 dz 


o o 


Tt/2 


B(K^) 


[(f W* + f*W ) cos 0 de] r, 
ss ss z=0 


7r/2 



• 

1 

~ r 


' « 

N2 

%■ 
— V 


3z 


W*i 



'F ^-1 

s 

+ QA 
s 


2H 

J 

L 3z 

2HJJ 


cos 6 d0 dz 


o o 


Thus, in the above equations the terms enclosed by angle brackets 

represent transfers of energy among the components K, A, K , and A while 

s s 

the terms G and represent diabatic heat sources and the terms B(K), 

B(A), and B(Kg) represent energy fluxes across the lower boundary. 

Summing (5.1) - (5.4) we see that the total energy K + A + K + A is 

s s 

conserved in the absence of diabatic heating and boundary fluxes. 



6. FINITE DIFFERENCE EQUATIONS 


6.1 The grid mesh 

All field variables are represented on a staggered grid in the 
meridional plane with grid points identified by the indices (j,k). Here, 
j increases southwards and k. increases upwards. To minimize truncation errors 
the grid points are staggered as shown in Fig. 1. The grid staggering is 
arranged in the horizontal so that U, V, X, 1” and W are defined at 
the meridional points given by y = (TTa/2) - (j - IjAy, j = 1,2,....J^ where 


Ay = \ 


TTa/ ( J -• 1) , global domain 
m 

, > , hemispheric domain J 

- i; 


m 


The variables 'f, U*, V', W and F^ are defined at the meridional points 

TTq 

y = - a - l/2)Ay, j = 1,2,....J^-1 . 

Thus, U, V, X, 'f* and W* are defined at the horizontal boundary points while 
Y, U' , V’, and W are defined a distance Ay/2 inside the boundaries. This 
form of staggering is natural for use with the horizontal boundary 
conditions (4.1). 

The vertical staggering is arranged so that U, V, T, U’ , V’, H*’ and 
are defined at the levels 

z = (k - l)Az, k = 1,2,....K^ 

where Az is the vertical grid increment; while the variables W, X, W’ and 
F^ are defined at the levels 

z = (k - l/2)Az, k = 1,2, * 
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6.2 The difference equations for the zonal mean 

For time differencing we choose a semi-implicit method in which the 
iner.tia-gravity terms (i.e., the Coriolis, pressure gradient, and adiabatic 
heating terms) are treated implicitly while the nonlinear advection terms and 
forcing terms are represented by centered differences. However, in order 
to prevent a weak time splitting of the solutions , which occurs due to the 
"leapfrog" scheme for the advection terms, a forward time step is used once 
every 48 steps. 

The time differencing scheme can be expressed efficiently if we define 
a time average as follows 

F = (6.1) 


Here F stands for any dependent variable, n is the time index given by 
t = nAt, n - 0,1,2 

where At is the time step, and B 2 » are defined as follows: 

(a) leapfrog step, = 1/4, $2 = 1/2, B 3 = 1/4 

(b) forward step, B^^ = 1 / 2 , = If 2 ^3 = 0 

For leapfrog steps the time difference can then be written as 


rwl” + l ^ 

8tJ 2 At (At/ 2 ) 


( 6 . 2 a) 


While for forward steps we have 


3t 




pU + 1 

At 


- F 


,n 


1? — T? 

(At/2) 


(6.2b) 
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In writing out the space differences it is convenient to use the 
following differencing and averaging operators: 


^l/aC ) = [( - ( )j+]^]/Ay (6.3a) 

^ +1/2 ^ ^ j +1^/^ (6.3b) 


Furthermore, to write the required vertical differences we let 


_F_ 
3z 2H 


= e 


•z/2H 9 z/2H. + -n in 

(Pe ' ) =. - F^a )Mz 


(6,4a) 


, + _ Az/4H j - _ -Az/4H , 

where e = e and e = e . Sxmilarly, we have 


- (F e - F e^)Az 
[9z 2 hJ '^^k+1 k '' 


(6.4b) 


where in each case the difference is centered at the k + 1/2 level.- 

Using the operators defined in (6.1) - (6.4) we can write finite 
difference approximations to (3.7) - (3.10) as follows: 


/s. /S 


U - (fAt/2)V = A 

(6.5) 

V+ (fAt/2)U + (At/2)6^ 

(6.6) 

(V - e + i/2(\ - 0 

(6.7) 

C¥ e"*" ^ e”) + ^ AtAz — ^ 

^\+l \ ^ 2 \ ^ 

(6.8) 


Here the terms involving the unknown variables have been collected on the 
left hand sides, and the source terms involving known quantities at time 
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levels n and n - 1 appear on the right hand sides . In writing out these 
eqxiations the subscripts j ,k have been omitted wherever no ambiguity would 
result. In the continuity equation cos 0 is required at both the V and 
W grid points. Thus, we define 

9^ = tt/ 2 - (j - 1/2) Ay/a (6.9a) 

0* = ir/2 - (j - l)Ay/a (6.9b) 

The source terms A. and B are defined as follows: 


A - ^ - ^ ((tj'^cos e*>j<v"cos 8*>.) 

+ 2A. Ls 6^ tWjV «>3 -l/2,k 

- + »VVa.k-l’} 


+ ^[F^ + Di(U° ^)3 


(6.10a) 


1 At z/2H 

& 

L 0 “ 3. 

cos 0* _ 

1 " 1 

cos 0* ■ 

+ ^ e 

4Ay. 

cos 0* 
1 

cos 0* T 
J - 


a. 


cos 0* 
1 

cos 0* . - -3' 
.1 + 1 

. At 

cos d*,. 
J + 1 

cos 0* 
1 

.. 


(6.10b) 


Here the coefficients and defined as follows: 

(a) leapfrog step, ]X^ - - 111 

(b) forward step, = 1, y^ = 0 
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In formulating B we have used a special form for the so-called "metric" 
term, tan 6/a, which is required to keep the difference equations energy 
conserving in adiabatic, frictionless flow. To derive this form we note that 


tan 0 




2 cos^ 0 dy 


XT (cos^ 0) 


This may be approximated as 


U. 

-JL 

4Ay 



fcos^ 6* - 

- cos^ 9*> 


rcos^ 6* - cos^ 9* , i'll 

”3-1 

3.-1 

X 

+ u. , . 
3+1 

3 .1+1 

cos 0* cos 0* , 

^ 3 3 - 1 ^ 

COS0* cos 0* , , 

J 3 + 1 


which easily reduces to the form given in (6.10b). 

In order to write the thermodynamic source term, R, in a compact form 
we define a density weighted "thickness" by letting 


„ n _ 


. T e"^ - e' 

[ 3,k + l j ,k 


( 6 . 11 ) 


We then have: 


i? = ^ [Az(F^ + Q) + ^)] 


^ z/2H j 1 . 

2 ^ \2 cos 0 j+1/2 






- (W? + w? )(S? e*^ + sj', 1 e )] 

j,k 3,k-l j,k 3,k-l 


( 6 . 12 ) 
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6.3 The difference equations for the eddies 

Using the above notation the equations may be written as follows : 


^ ^ -im At 

% - (£At/2)V^ = 




(6.13) 


\ + (fAt/2)S^ = f s.^y^ 


(6.14) 


, N^AtAz r. 

('P 1 . 1 e - 'F , e ) + — z W = E 

s,k+l s,k 2 s s 


(6.1-5) 


CW, 




Where here = s/(a cos 0) 


. n 




')+% s coa e “3+‘/2 




«a ^ Vs"" - f VVa 


(6.17) 
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s 

u. 

J 

fcos 0* 
3 

cos 0.-V 
3 

Ay 

cos 0. 

cos 0* 


^ 3 

3-^ 


+ U 


j+1 


cos 0. cos 0 *.t 1 

J 1±1 


cos 0*,, cos 0. 
3+1 j 


jj 




(6.18) 
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(Note that if U is symmetric about the equator then A and B have the same 

s s 

symmetry as U and V , respectively.) 

^ s s 

In formulating B ^ we have used a special form for the metric term, 

-2UU tan0/a, which is required to keep the difference equations energy 
s 

conserving in adiabatic frictionless flow. To derive this form note that 


H 

2' U U tan0/a = (cos^ 0) 

s cos-^ 0 dy 


This may be approximated as 


U r /-COS^ 0* - COS^ O.’v 

^ u J 1 

Ay i cos 0* cos 0 . 

^ ^ ^ j 2 ^ 


rcos^ 6. - cos^ 0*. 

j 3 , 

i+1 I cos 6. cos 

^ 13+1 


which easily reduces to the form given in (6.18). 

The term 2?^ in (6.15) can be written in a fairly compact form if we 
first define a density weighted "thickness" 


S . = QH . , r 1 - 'i' . 1 e ) 

s,l,fc s,j,k+l s,3,k 


We then can write 


R = y S “ Q + ^ D (£T ’^ ^ ) 

s *^1 s,j,k '^2^3, j,k 2 2 2 s 


At z/2H n — n ^ „n 


w 

Az/2H y^n\ -Az/2H 

+ 2 Az p ^j-V 2 .k+i = = J 


). 1/ v+1 + )• 

'■ S '.j- / 2 ,k+l ' « 't- 


■'i-Va.kl 5 , jn 

2 J ‘j-Vz ^k 



7. SOLUTION METHOD 


7.1 The zonal mean equations 

The system (6.5) - (6.8) is a set of simultaneous equations for the 

/S /N /N /N 

unknowns U, V, and W. To solve this set we first eliminate U between 
(6.5) and (6.6) to obtain 


v = y. 




(7.1) 


where y. = (1 + f^At^/4)’’^, 

y*v 

We next substitute from (7.1) and (6.8) into (6.7) to eliminate W 

/N 

and V. The result is an elliptic difference equation in H*: 


V- v+i - 1 f, ^ 

k j ,k+l k-1 k jjk 


+ r, , 4'. . T + A. 'F. , . + B. 'F. , + C. = D. , 

k-1 j,k-l j J-l»k J J,k j 3+1, k j,k 


(7.2) 


Here we have let N^(z) = H^/r(z) where N^ = constant, and then expressed 

r(z) (the vertical variation of static stability) at z = (k - 1/2) Az as F 

The coefficients. A., B., C., D. are then defined by 

3 3 3 3 ^ 

N^Az^At^ 

a (Y- cos 0*) 

3 4Ay cos 0^ 'j 2 


k' 


N^Az^At^ 

S " AAy"" cos e. ^"'^ 3+1 ®j+l^ 


B. E - A. - C. 
3 3 3 
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k: 3>k 


r, T i? . , T 

k-l 3,k-l 


+ 

e 


N^AtAz^ _ _ 

+ f^STeT VV2 ® - 

The elliptic system (7.2) may be solved for ¥. , using the NCAR 

J > 

subroutine BLKTRI (Swarztrauber and Sweet, 1975). In this case the solution 
is carried out for the grid points in the range 




The lateral boundary condition (4.1) is incorporated by setting 


= 0 and “ 0 

m 

✓N 

The lower boundary condition is incorporated by setting T. - equal to the 

J » -*- 

value obtained from integrating (4.2c). The upper boundary condition 
(2.16c) in finite difference form requires 




e 


(7.3) 


Once has been obtained by inversion of (7.2) it is a simple matter 
to compute the remaining fields. If, however, one attempts to compute V 
from (7.1) the results are rather unsatisfactory due to large truncation 
errors. This problem arises due to the fact that the first and third terms 
on the right side are generally two orders of magnitude greater than V so 
that V Is obtained as a small residual of two large hut opposite terms. To 
avoid this problem it is useful to utilize the meridional streamf unction 
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defined by 



(3.11). 

1 

cos 6 


In finite difference form we have 



X 


(7.4a) 



1 

Az cos e* 



3 



(7.4b) 


which identically satisfies the finite difference form of the continuity 


equation (6.7). 


Substituting from (7.4a) into (6.8) and noting that = 0 we can 

/S 

solve for X. , : 

J 


• = V -2Ay cos 9 

•j+l,k j,k H^^zAt ^ 




e )] 


(7.5) 


We next use (7.4b) to solve for V and finally use (6.5) to obtain U. The 
final step of the solution is then to use the definition (6.1) to obtain 
all fields at time n + 1. For example. 


= (U - ^2^^ - 

and similarly for , and . 


(7.6) 


7.2 The eddy equations 

The system (6.13) - (6.16) is a set of simultaneous equations for the 




unknowns U^, V, W^ which is exactly analogous to the mean flow set 

/N 

discussed above. To solve this set we first solve for U and V in terms of 

s s s 


> 
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using (6.13) and (6.14): 


U = Y. 

s ' j 


r-im At ^ ^ 


f2At 


6 . 


^ + A + ^ B 


4 ^j+V, ^ 2 "s 


(7.7) 


V = Y- 
s 2 


At 


2 , ^ 


fAt 


im f ^ ('5',)-. 1/ - 5. . 1/ '1' + B 

, s 4 s j+ /2 2 j+ /2 s s 2 s 


(7.8) 


where y^ = (1 + f^AtV4)“\ 

Combining (6,15), (6.16), (7.7), and (7.8) we get a single equation for 


tJ = (r ^ j, -Az/2H^ J + T ¥ ■ 

k Sj k-1 k ^ s,j,k k-1 s,j,k-l 


+ D .¥ ,,,+E .¥ .,+F Y .^-,,=T (7.9) 

s,j s,j-l,k s,j s,j,k s,j+l,k s,j,k 

where T, is as defined below (7.2), and the coefficients D , E , F are 

K. S S S 


N^Az^At^ r-Y. T cos e. ^ s. ^ ^j-1 ®j-l-i 

D =_2 -lul ^ 

s,j 4 cos 6* Ay2 4 

2 


N^Az^At^ 


F . 

s,J 


_ o 


4 cos 0* 
2 


Y. cos 6. Y. cos 0,-r 


Av^ 


J. _ s j 


iL 


E . 


Az^ At^ 



4 cos 0* 

2 


m 


Jzi. 


fY . , cos 
■1-1 

®.1-i "'.i 

. cos 0 , 

_ 1. . 1 


Ay2 


i-i 

+ m cos 
s . 

1 

0. Y- 

2 2 

4 


i — :: — (y. T m cos 0. - f. , ~ y.^i cos 9, f.) 
Ay 'j-1 j-1 j-1 'j s^ 2 2 
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and the source term is 


^s,j,k ^^k^s,j,k® ^k-1 -^s,j,k ® ^ 


At Az^ 

[(Y- ^ cos 6. T q . -, - Y. cos 6. q .)/Ay 
2 cos 8* j-1 'j j 


+ ^ (m . - Y . - cos 0 . - p , - + m . y . cos 9 . p . ) ] 
2 s,3-l J-1 J-1 s,j-l s,j 'j 3 ‘^s.j 


where 


, , fAt „ n, fAt - 

p .= A + -:r- B , q . ~ B t- A . 

s ,3 s 2s’ ^s,j s 2 s ,3 


The elliptic system (7.9) may be solved for 'f . , using the NCAR subroutine 

s > J »k 

CBLKTRI. 

For global or hemispheric antisynmetric modes the solution is carried 
out for grid points 


2<j<J^-l; 2<k<K^-l 


However, for symmetric hemispheric calculations the solution includes the 
point 3 = J^. 

In the global or antisymmetric hemispheric case we thus require 


^ T V = 0 

s,l,k s.J ,k 

’ ’ ’ m 


(7.10) 


while for the symmetric hemispheric case we must have ^ , = 0 and 

s , 1 , k 


'F 

s,J T,k 
m-1 


= 'F 


(7.11) 
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while 


Condition (7.10) is incorporated by letting D -j = 0 and F , =0 whi 

s , 2 s , Oni-l 

the condition (7.11) requires replacing D ^ as defined above by D ^ + F 


In all cases the upper boundary condition is Y , = 0 and the 

S 

lower boundary condition is a specified forcing 


. -I = gh(y,t) 

» s J * -L 


(7.12) 


Once ¥ , is obtained from (7,9) we compute W . , from (6.15) 


^s,j,k Atir^" Az ^■®s,j,k ^ '^s,j,k ® 


(7.13) 


We then use (7.7) and (7.8) to solve for U and V . Finally, these results 

s s 

are used to obtain ^ formula analogous to (7.6). 

A similar treatment of proved unstable. Therefore in computing fluxes and 


vertical advection terms W is used in place of W 

S 5 


n+i 


7.3 The eddy flux terms 

The eddy momentum flux convergence (3.11) and the eddy heat flux 
convergence (3.12) must be written in finite differences so that the energy 
integrals of the flow remain satisfied. It turns out that energetically 


consistent forms are: 
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F = - e 
m 


z/2H 


j <s. 1 , t(u V* + U*V ) cos^ 9 ] 

[cos-^ O'? 3 - 72 s s s s 


K^sj.k °so.k+l 


- <U . T e'*' + U . . , e ). 1 , W* . , , 
' s,j,k s,j,k-l /j- /2 s,3,k-l 


^^s,j,k ® ^s,j,k+l ^ )j-V2 ^s,j,k 
- <"13, k + "2,3.k-l "')3-Va «.,3,k-l'} 


(7.14) 


and, 


(z+Az/2)/H 


‘ 2 cos 6. "j+Y 2 I 


6,,i. -^cos 0* [(<V ). 1, + (V ). 1/ , ) S* . , 

J+/2 I 3 ' s'j-72,k+l ^ s'j-V 2 ,k s,j,l 


^ «''^>3-Va,kkl + <''?)3-'A.k> ^,3.k'} 


zIlE 

2Az *-^^s^s^J4-V2,k-4-l “ ^^s^s)j+V 2 »k-l 

^ VV2,k-M - 


(7,15) 



8. INTEGRAL CONSTRAINTS AND SUBGRID SCALE DIFFUSION 


8-1 Integral constraints for the zonal mean equations 

The basic equations of the model (3.7) - (3.10) satisfy certain integral 
constraints which also must be satisfied by the finite difference equations 
if satisfactory long term integrations are to be obtained. It is easily 
verified that when the forcing terms F^, and Q are omitted, and sub grid 
scale diffusion is neglected, the rate of change of zonal mean kinetic 
plus available potential energy is equal to the energy flux through the 
lower boundary: 

A (P + K) - j h g + ^ -- 

A 


W 


z=0 


dA 


( 8 . 1 ) 


where 


P 


1 

3'1' ¥ ' 

2 

J 

3z 2H 


T 


/N^dx 


(8.2a) 


K = 



T 


(u2 + v2)dx 


and 


(8.2b) 


dA = a^ cos 6 d6 dA 
dx = a^ cos 0 d6 dA dz 


Another important constraint is the conservation of relative angular 
momentum. If we multiply (3.7) by exp (-z/2H) cos 0 and integrate the result 
over the entire domain we find that relative angular momentum is conserved 
except for the flux of angular momentxun through the lower boundary: 
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_d_ 

dt 


-z/2H - 

e If cos 0 dT = 


[U W cos 6] ^ dA + 

z=0 


fX(z=0)dA 


(8.3) 


In deriving (8.3) we have neglected eddy momentuin fluxes through the lower 
boundary. It is important to note in connection with (8.3) that horizontal 
diffusion can not change relative angular momentum so that we require 


j D^(TJ) cos 0 dA = 0 (8.4) 

A 

which constrains the possible forms for the operatore D^( ). 

Also, horizontal diffusion can not change the horizontally averaged 
temperature (thickness) on a horizontal surface. Thus from (3.10): 



'3'!' , 


J 2 

3z 

2H 


A 


dA = 0 


(8.5) 


The constraints (8.1), (8.3), and (8.5) must also be satisfied by our 

2 

system of finite difference equations if satisfactory results are to be 
obtained. In finite difference form the integrals are replaced by sums 
over the grid points : 


r 


J 


( )dx ^ 


L 




( )2TTa cos 6 Ay Az 


(8.6a) 


( )dA ^ ( )27Ta cos 0 Ay 
A j 


Except for the effects of time truncation errors. 


(8.6b) 
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The -value of 0^ used in (8.6a) or (8.6b) is given by either (6.9a) or 
(6.9h) depending on the location of the dependent variable; e.g., in (8.5) 
we use (6.9a) and in (8.3) we used (6.9b). 

Next, multiplying (6.5) by e cos^ 0*(2iTaAyAz) we find after summing 

over all grid points (2 < j < - 1, 2 < k < - 1) , and using 6.2a): 

^ I U cos^ 9* = X ^1.1 

j A j 


e ) <W cos cos 0*] (8.7) 


which is consistent with the differential form for angular momentvun conservation, 

(8.3). (We have here assumed that (8.4) holds for the finite difference 

form of subgrid scale diffusion.) 

The finite difference analogy to the energy integral (8.1) may be 

« 

obtained by multiplying (6.5) by U cos 0*, (6.6) by V cos 0*, and (6,8) by 
(cos 6 )/N^Az^) ('l'^_j_^ 'e"*" - e ) then adding the three resulting equations 
together and summing over all gridpoints. Using (6.2) to express the time 
derivatives in differential form we can then write 


(U?, + V?, )cos 0* 
dk ik 1 


IAz^N^ 


j=2,J -1 
•J ’ m 

k=2,%-l 


3 '^m 

k=l,N 


E r- - - 9 U 

W. _ 'F. e cos 0. d li— (y ^os 0), i, . e 

3 2 ' ^ 3 - / 2,1 


3 “1 » 3^m 


4N^ Az^ 


e 9 e -1 
Lit 


(8.8) 
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where = Q = 0 and friction terms have all been neglected. Clearly, 

(8.9) is a reasonable analogue to the differential relationship (8.1). 


8.2 Sub grid scale diffusion for the mean flow equations 

In order to suppress nonlinear instability it is necessary to smooth 
all fields in the meridional direction. In order to prevent this smoothing 
from damping the large scale motions we have chosen to use a fourth order 
linear diffusion operator. In applying the diffusion in the zonal momentum 
and thermodynamic energy equations we must recall that both relative angular 
momenttim and horizontal average temperature must be conserved [see (8.4) and 
(8.5)]. In addition the diffusion terms should make negative definite 
contributions to the energy equation. 

In order to satisfy both these requirements it turns out that in the 
zonal momentum equation relative angular velocity should be diffused. Thus, 




K t" 
cos^ 0 0y2 


cos 0 

V / 


(8.9) 


This automatically satisfies (8.4) provided that (8 ^/3y ®) (U/cos 0) = 0 at 
the meridional boundaries. If we multiply (8.9) by U cos 6 and integrate 
the result in y we obtain 


• _ _ 

D cos 0 Dj^(U)dA = - K 

A 


[3" 

r u 

by" 

cos 0j 


dy dX 


which is negative definite. Thus, the diffusion term (8.9) acts as an 
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energy sink. In finite difference form we write 


D^(U) 


- K 

rr u ^ 

— A 

f ““ ^ 

U 

cos 2 9* Ay'* 

[cos 0* 

3-2 

cos 0* 

V. J 


u 

A 

f > 

u 

4* 

full 

cos 9* 

. H ; 

J 

cos 0*j 

j+1 

cos 6*Jj+2 


( 8 . 10 ) 


In order that this finite difference form satisfy the difference analogue 

of (8.4) i.e., cos^ 6* I>j^(U) = 0 the formula must be modified at the 
3 

points adjacent to the boundaries. Thus 


[D^(U)] 


j=2 


- K 

'■ 2U ' 


3U ' 


U 

1 1 

cos^ 9* Ay** 

cos 9*^ 

3=2 

cos 0*^ 

- -5 ^ 

3=3 

cos 9* 

k J 

It 


(8.11a) 




- K 


0* Ay‘ 


U 


cos 0 


3=2 


+ 6 


U 


cos 0* 


3=3 


- 4 


U 


cos 0* 


3=4 


4" 


U 


cos 9* 

J 

3=5 _ 


(8.11b) 


with analogous expressions for j = - 1 and j = - 2. Again using the 

notation of (6.11) we can write the finite difference diffusion term in the 
thermodynamic energy equation as follows: 


(S) — O A 4 O Tw , ^ 

2 cos 0 Ay^ j-2 j-1 


IS. , - 4S_ + 6S, - 4S._^ + S.^^] 


( 3 . 12 ) 


Again the points adjacent to the boundaries require special treatment: 


[D.,(S)],_, = 


- K 


2"‘""-'j=l cos 0 Ay'* ^^“j=l -^“j=2 3=3 




(8.13) 
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[D_ (S) ] . - = ~ :r . [-3S . , + 6S . „ - 4S . . + S . , ] 

2 3=2 COS e Ay“* j=l 3=2 3=3 3=4'' 


(8.14) 


K 


^^2^®^^j=Vl cos 0 Ay“ 


(8.15) 


Finally we write 


^ 2 ^ = cos'el Ay^ ^ - «3+l + ''3+23 


(8.16) 


with the special cases 


[D„ (V), 


- K 


2 ^‘'^ 3=2 “ cos 0 * Ay'* '"’' 3=2 ^''j=3 " ''j=4 




(8.17) 


and an analogous form for 3 = - 1 . 


8. 3 Subgrid scale diffusion for the eddy equations 

To filter out small scale noise so as to suppress nonlinear instability 

the eddy equations include fourth order linear diffusion terms similar to 

those discussed in Section 8.2 for the zonal mean flow. D„(U ) and D (V ) 

^ s ^ s 

have the same form as 1 * 2 ( 3 ) given in (8.12), while ^ 2 ( 1 ?^) has the form of 

^ 2 ^^^ given in (8.16). These forms must, however, be modified next to the 

boundaries to insure that diffusion does not change the meridional average 

of any field. For a global domain the modification to D^(5 ) is identical 

to that given in (8.17) for D 2 (V). 

For the U and V field the situation is more complicated since the 
s s 

boundary conditions are different in the s = 1 and s > 1 cases . 
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For the case s = 1, D (U ) and (V ) are computed using formulas 

Z s ^ s 

analogous to (8.12), (8.13), and (8.14). For global integrations formulas 
similar to (8.13) and (8.14) are applied at j = '7^-1 and j = - 2. 

For the case s > 1 the polar boundary condition requires that 


h,(U = 


K 


[^U - 3U + U 


2" s"j=l cos 6 Ay** ^ s,j=l s,j=2 s,j=3 


(8.18) 


^2(^s>j=2 = 


K 


[-5U . , + 6U . „ - 4U . „ + U . (8,19) 


S»J=1 s,3-2 s,j=3 s,3=4 


with similar expressions for . 

In the case of hemispheric integrations the diffusion operators at 

i = J - 1 and i = J - 2 are modified as follows : 
m m 

For antisymmetry conditions on ¥ the form of D„ (5 ) is the same as in 

^ s 

the global case; however, since TJ = - U , and V = + V ^ ^ 

S,J^ SjJjjj ± -L 

we use a diffusion form analogous to (8.18) and (8.19) for D (U ) ^ - and 

L s Jin“i 

D (U ) n while for D (V ) ^ ^ and D„ (V ) ^ „ we use forms analogous to 

® ‘^m ^ 2 s Jjji i 2 s 

(8.14) and (8.15). 

For symmetric conditions on '1' the form of must be modified as 

follows : 


D.C-S ). •, = T +75, 1 - 9 +^T 

2 s J -1 cos 0 Ay^ s,J S,J -1 J -2 J -3 

m ’mm m m 


- K 


^2^^s^J cos G Ay 
m 


IT E35j - 45 
m 


J -1 - 7 } 

m m 


( 8 . 20 ) 

( 8 . 21 ) 


and since U _ = U ^ ^ and V ^ = -V _ , we use forms similar to 

s»Jm ®’-7m ^,1^,-1 

(8.14) and (8.15) for D (U ) . and D (U ) and forms analogous to (8.18) 

2 s ^ ^ 

and (8.19) for D_ (V ) ^ ^ and D (V ) ,, 

2 s Jm-1 2 s Jjq-*2 



9. DIABATIC HEATING COMPUTATION 


9.1 Infrared heating 

This study has utilized Dickinson’s (1973) parameterization of infrared 
cooling consisting of the sum of the cooling for a reference temperature 
T^ and a Newtonian cooling approximation for the departures from that profile. 
Thus the net heating teinns take the following forms : for the eddies , 

Q = - aT 
s s 

while for the mean flow 

Q = - (Q^ + oT) 

where is the diabatic heating due to the absorption of solar radiation by 
ozone. is the net infrared cooling at each level for the reference 

temperature profile, and a is the Newtonian cooling coefficient. and T 
are functions of altitude and latitude while a, and T^ depend on 
altitude alone. 

The values of the Newtonian cooling coefficients have been calculated 
for levels between 30 and 80 km by Dickinson (1973). Below 30 km Trenberth’s 
(1973) values are adopted. Although the accuracy of the Newtonian cooling 
representation breaks down above .about 70 km, it shall be retained at this 
time for lack of a better representation. Following Schoeberl and Strobel 
(1978), the value of a between 80 and 96 km was taken to be the .CO^ cooling 
rate in the fundamental band at 15u (see Fig. 2). 

Dickinson’s (1973) careful computations of a and were made for 
atmospheric temperature profiles that differ little from the reference 
temperature profile. Because the actual temperatures may vary considerably 
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from this reference profile, especially in the winter polar region, alternative 
values of are here computed in the following manner. 

At a given level the globally averaged diabatic heating Q is given by 

Q = Qg “ (Q^ - aT) 

where (~) designates a horizontal average on the sphere. Since the observed 
globally averaged ten^jerature profile is fairly well known, we choose 

so that global radiative equilibrium (Q = 0) is achieved when the globally 
averaged temperature profile T is equal to T^. Therefore 



9.2 Solar Heating 

Below 96 km ozone is the only significant absorber of solar radiation. 
The parameterizations of Lacis and Hansen (1974) are used to cot5)ute the 
solar heating term Q^. The diurnally averaged solar heating is calculated by 
fixing the sun angle at its average value between sunrise and sunset 
(approximation 1 of Cogley and Borucki, 1976). The sun angle may remain 
fixed for the duration of a given rmi, or may be varied according to the 
seasonal cycle depending on the objectives of the particular run. 



10. A TEST APPLICATION OF THE MODEL 


In order to demonstrate the capabilities of the model we have computed 
the zonal mean annual cycle for the stratosphere and mesosphere for conditions 
of zonal mean forcing only. In this experiment the eddy forcing was set to 
zero at the lower boundary. The mean zonal winds at the lower boundary 
(16 km) were specified to vary over the annual cycle according to the 
observations of Labitzke (1972) for the northern hemisphere and Taljaard et al . 
(1969) for the southern hemisphere. The diabatic heating was also specified 
to vary on the annual cycle by including seasonal variations in the solar 
zenith angle and sun-earth distance. 


10.1 Rayleigh friction parameterization 

In order to produce a realistic mean wind profile it proved necessary 
to specify strong damping in the mean momentum equations above 70 km. In 
the atmosphere the mechanical damping of the mean wind near the mesopause 
is probably due to the breaking of gravity waves and tides. For the present 
model this effect is parameterized in the simplest possible form by using 
a height dependent Rayleigh friction coefficient 


K. 


R 


K 

o 


+ 1C, 


1. + tanh 



where ic^ - 1/80 days, = 1/4 days and z is in kilometers. This profile is 
shown in Fig. 2. 

The biharmonic horizontal diffusion coefficient is given the value 
K/Ay"* = 10“® which is the minimum necessary to suppress nonlinear 
computational instability when At = 1 hr. 
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10.2 The zonal mean annual cycle 

Figs. S-'S show the zonal mean wind, mean meridional wind, and vertical 
velocity profiles for .southern hemisphere winter solstice and spring equinox 
conditions computed using the above described parameters and a grid resolution 
of 10° latitude and 5 km height. During the solstice season there is a 
thermally direct mean meridional circulation with rising motion in the summer 
hemisphere and sinking in the winter hemisphere. At the equinox, on the 
other hand there is a two cell meridional circulation with rising in the 
equatorial zone and sinking near both poles. Zonal mean winds computed in 
both seasons are quite realistic. This example shows that a zonal mean 
model is capable of simulating many important features of the general 
circulation in the middle atmosphere. Further details of this annual 
cycle simulation are reported in Holton and Wehrbein (1979) . 
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Figure 1; A portion. 

of the 

grid mesh in the meridional plane showing 

the arrangement of variables on the staggered grid. 



Figure 2; Vertical profiles of the Newtonian cooling coefficxent 
(solid line) and Rayleigh friction coefficient (dashed 
line) in units of d*”^. 
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Figure 3: Computed mean zonal winds (m ) for the Southern 

Hemisphere winter solstice. 
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Figure 4: Computed mean zonal winds (m s“^) for the Southern 


Hemisphere vernal equinox. 
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Figure 5: Computed mean meridional wind (m ) for the 

Southern Hemisphere winter solstice. 
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Figure 6: Computed mean meridional wind (m s~^) for the 

Southern Hemisphere vernal equinox. 
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Figure 7 : Computed mean vertical velocity (mm s“^ ) for the 

Southern Hemisphere winter solstice. 
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Figure 8; Computed mean vertical velocity (mm for the 

Southern Hemisphere vernal equinox. 




APPENDIX 


FORTRAN CODE FOR THE SEMI-SPECTRAL MODEL 

The present version of the model computes the interaction of a single 
wave mode with the zonal mean flow. The program consists of the main 
program, PROGRAM WA^7E2, in which the fields are initialized and the calling 
sequence for the various subroutines is established. The main dynamical 
computations are carried out in SUBROUTINE ASTREAM (mean flow equations) 
and SUBROUTINE EDDY (wave equations) . The radiative heating calculations 
are carried out in SUBROUTINE' HEAT , SUBROUTINE RADEQU, FUNCTION BELT, and 
FUNCTION OZUV. The output fields are created in SUBROUTINE AOUT. All 
of the above routines are given in the following program listing. In 
addition the program requires the NCAR library routines SUBROUTINE BLKTRI 
and SUBROUTINE CBLKTRl. 

In order to facilitate reading of the code a dictionary of the principal 


FORTRAN sysmbols is provided below. In addition to defining the symbols, 
we have Indicated the location in the program where each symbol first appears. 



Dictionary of FORTRAN Symbols 


FORTRAN Symbol | 

Definition 

Location 

1 

1 

A 

a, Newtonian cooling coefficient 

C 24 

AIRDEN 

air density 

C 36 

ALA(J> 

At/Ay (2Q sin 0) 

A 71 

ALB(J) 

At(2S^ sin + 

A 72 

AM(J) 

A. coefficient in (7.2) 
3 

A 103 

AN(K-l) 

/\ 

Coefficient of ^ in (7.2) 

A 175 

ANG(J) 

s U/(a cos 6*) 

G 74 

AR(K) 

Coefficient of 1st term in (7.9) 

A 180 

AS 

A (6.10a) 

F 106 

AZEN 

Magnification factor 

C 38 

BLKTRI ' 

Subroutine to solve (7.2) 

F 139 

BM(J) 

B. Coefficient in (7.2) 
3 

A 106 

BN(K-l) 

Coefficient of 2nd term in (7.2) 

A 177 

BR(K) 

Coefficient of 2nd term in (7.9) 

A 182 

BS 

1 

B (6.10b) 

F 113 

BVF 

(buoyancy frequency) 

A 29 

B1(J) 

E . coefficient in (7.9) 
Sj3 

A 188 

CC02 

Infrared heating 

C 25 

CDUM(J,K) 

Dummy Array 

G 42 

CG(J,K) 

(unused) 

A 42 

CHI(J,K) 

R, ^ (6.12) 

F 61 

Cl 

A 52 

CLOUD 

fractional cloudiness 

C 41 




FORTRAN Symbol 

Definition 

Location 

CM(J) 

C. coefficient in (7.2) 
1 

1 • 

A. 105 

CN(K-l) 

Coefficient of , in (7.2) 

k+1 

A 176 

CNA(J) 

N^Az^At/(2Ay cos 0) 

A 77 

CNB(J) 

At/(2Ay cos 0) 

A 68 

COR(J) 

2fi sin 6 

A 65 

CORIOL(J) 

1 2J2 sin 0* 

A 66 

COZONE 

(unused) 

C 8 

CR(K) 

Coefficient of 3rd term in (7.9) 

A 180 

CS(J) 

cos 0 

A 63 

CSA(J) 

cos 0* 

A 64 

CTS(J,K) 

(unused) 

A 42 

CUB(J) 

y cos 0* i 

A 78 

CURTIS (J,K) 

(unused) 

A 37 

C1(J) 

F . coefficient in (7.9) 
j 

A 187 

DCA(K) 

a, At/2 
k 

A 227 

DCOS(J) 

1. /(a cos 6*) 

A 87 

DELAY 

forcing switch-on time delay 

A 229 

DEL(J) 

N^AtAz^/cos 0* 

A 85 

BELT 

At 

A 30 

DELTA ■ 

solar declination 

B 24 

DENS (K) 

exp(z/ 2H) 

A 94 

DKY 

K/Ay"* 

A 228 

DR(J) 

B - JfAt/2 

F 130 

DT 

Lull 

A 31 





FORTRAN Symbol 

Definition j 

Location 

DTODY 

At/ (4Ay) 

A 57 

DIM(J) 

dtumny 

I 18 

DY 

Ay 

A 51 

DY2 

(Ay)^ 

A 54 

DZ 

Az 

A 29 

DZ2 

Az^ 

A 53 

DZN 

N^Az^At^/ (4Ay^ cos 0) 

A 103 

EC 

Orbital eccentricity of the earth 

B 21 

EMI 

exp(-Az/4H) 

A 59 

EPL 

exp(+ z/4H) 

A 58 

FDAY 

fraction of day that sun shines 

B 49 

EM(J,K) 


H 38 

FREQ 

Q 

A 55 

FT(J,K) 

(7.15) 

H 22 

GAM(J) 

\ 

A 67 

GAMB(J) 

— ^ 

A 76 

GBV(K) 


A 95 

GMEP(J) 

vertical advection of T 

F 76 

GMl(J) 

.51 s At/(a cos 6) 

A 89 

GTIME 

growth time for forcing 

A 231 

ICLOUD 

altitude layer of clouds 

C 43 

ICT 

Index for forward difference 

A 33 

lEND 

Total time steps for run 

A 45 

IFD 

Frequency of forward steps 

A 32 

IFLG 

IFLG = 0 initializes BLKTRI 

A 29 







FORTRAN Symbol 

Definition 

Location 

IMAT 

Index for output 

A 111 

INIT 

Input Flag 0 for TAPE 1 input on 

continuation runs 

A 45 

IPHAS ( J) 

Dummy 

I 89 

IPRINT 

Frequency of output 

A 45 

TRAD 

Index for calls to HEAT 

A 99 

IRCT 

Frequency of calls to HEAT 

A 100 

ITIME 

Index for time step n 

A 110 

JM 

J 

m 

A 29 

JML 

J - 1 
m 

i A 47 

KAP(K) 

a^, Newtonian cooling 

A 

97 

KN 


; A 29 

KNL 


A 48 

KZ 

Diffusion Coefficient 

A 34 

M 

J - 2 
m 

A 49 

N 


A 50 

NGC(J) 


A 81 

NUl(J) 

.5i s(y cos 0)/(a cos 6) 

A 83 

OZUV 

Solar energy absorbed by ozone 

D 

1 

P(J) 

P 

s,j 

G 140 

PB(J,K) 


A 159 

PBA(J,K) 

Dummy array 

A 138 

PBO(J,K) 


A 160 

PERH 

Date of herhelion after Vernal Equinox' 

B 

18 

PHBI(J) 

Amplitude of boundary forcing for wave 

A .69 






FORTRAN Symbol 

Definition 

Location 

PHI 

i 

Latitude 

B 30 

PI 

ir 

A 46 

PLl(J) 

Y (k-1) ! 

S 

A 234 

PRS(J,K) 

Contains Fourier Coefficients of Dg(y,t) 
on input 

A 39 

P1-(J,K) 

^n + 1 
s 

A 243 

PU(J,K) 

Dummy array 

A 239 

P10(J,K) 

" ,}f n 
s 

A 243 

Q(J) 

s,3 

G 141 

QA(K) 

Globally averaged net radiative beating 

A 259 

QB(J,K) 

Q, diabatic heating 

B 93 

QDOT( ) 

Net radiative heating function 

B 89 

Q0(J,K), QOG(J) 

Ozone density/Lochschmidt’s number 

A 40 

QOZS(J,K), QOSZG(J) 

Ozone column abundance 

A 41 

QR 

Radiative cooling of reference atmosphere 

C 23 

QRS (K) 

Globally averaged solar heating 

A 115 

QS(J,K) 

Q , diabatic heating in wave 
s 

A 195 

R(J,K) 

T (on input) 

1 ® > J 

G 143 

RAB 

Effective albedo of lower atmosphere 

C 49 

RAD 

a (radius of earth) 

A 29 

RAYF(K) 

K_j Rayleigh friction 

A 225 

EEl 

Albedo of reflecting region 

C 51 

RDB 

Spherical albedo of reflecting region 

C 50 

RED(6,K) 

Newtonian cooling parameters 

A 43 

RG 

Ground reflectivity 

C 42 





FORTRAN Symbol 

Definition 

Location 

RHO 

Distance to sun in A.U. 

B 22 

RHZ 

H/287.1 

A 56 

RR(J,K) 

D (on input) 

J 5 k; 

F 135 

S 

Planetary Wavenumber 

A 26 

SA 

6* 

A 62 

SB 

e 

A 61 

SH 

H 

A 29 

S03 

1 

Ozone UV heating 

C 

63 

STAB(J) 

N^At^Az^/cos 0* 

A 86 

STRDAY 

Starting day (vernal equinox = 80) 

1 

A 36 

SU(J) 

sin 0 At^y/Ay 

A 80 

SV(J) 

YAt/Ay 

A 79 

SZT(J,K) 

Ozone UV heating field 

B 

67 

T(J,K), TG(J) 

Standard Atmosphere temperatures in each 
zone 

A 38 

TAU 

Optical depth of clouds in visible 

C 40 

TB(J,K) 

T, deviation of zonally averaged 
temperature from global mean 

1 

B 

76 

TEMP 

Dummy 

A 147 

TIME 

t = nAt 

A 113 

THETA 

Polar angle 

B 59 

TN(J) 

(4Ay)“^(cos 0*. -/cos 6* - cos 0*/cos 0* -) 
3-1 J 3 3-1 

A 91 

TNA(J) 

(cos 9*/cos 0, - cos 0./cos 0*)/Ay 
3 3 3 3 

A 82 

TNB (J) 

(cos 0./cos 6*,- - cos 0* -/cos 0,)/Ay 

j j+1 j+1 3 

A 90 

TR 

Reference temperature profile 

C 22 

TSTAR 

Local time of sunset 

B 40 





FORTRAN Symbol 

Definition 

Location 

TUl(J) 

i Y s At/(2a cos 6) 

A 74 

TVl(J) 

i Y s At^ (2fi sin 0)/2a cos 9) 

A 75 

TZO(K) 

Global radiative equilibrium temperature 

B 61 

T1 


F 25 

T2 

^2 

F 26 

T3 

^2^^1 

F 27 

T4 

33 / 6 ^ 

F 28 

T5 

l/3j^ 

F 29 

UB(J,K) 


A 127 

UBO(J,K) 

u’" 

A 128 

UT 

ozone path 

C 57 

U1(J,K) 

s 

A 243 

U10(J,K) 

s 

A 244 

VB(J,K) 

— n + 1 

A 250 

VBO(J,K) 


A 250 

V1(J,K) 

v”+‘ 

S 

A 244 

V10(J,K) 

V 

s 

A 244 

■WB(J,K) 


A 250 

WBO(J,K) 


A 250 

WA(I) 

Work Array 

A 249 

NR0(I) 

Work Array in BLKTRI 

A 244 

¥1(J,K) 

w"+' 

S 

A 244 

XBA(J,K) 

X 

A 250 

XMAl(J) 

s/(a cos 6*) 

A 88 





PORTKAN Symbol 

Definition 

Location 

XMl(J) 

s/ (a cos 6) 

A 73 

Z(K) 

z 

A 93 

ZEN 

Average value of cos (solar zenith angle) 

B 37 

ZT 

z + Az/2 

I 31 








c 


c 


c 


S$SS$$$$$(£$$ WAVE ZONAL FLOW INTERACTIO.N $$$SJSSSJtS4S$J 
PROGRAM WAVE2 ( INPUT» 0UTPUT>TAPE5-IHPUT#TAP66-DUTPUT»TAPei»TAPE2, 
$TAPE3) 

COMPLEX U1119>17)»U1D(19»171» Via9»17)»VlOa9»17)»Pl(19#171>P10(lS 
t, i7),Wl(19>16 >,W10C19i.l6),PlA{19,17)» Bia7>>R{17>15>>P(19)>QC19), 
JNU1(19)/CI>GK1(19)»TU1 (19)^ TVi(19 »» CDUM (19» 16 ) > PLl { 19 ) # QS (19/ 16 ) 
COMPLEX Y2/PHBKC19) 

PEAL UB(19/17J> U80(19/17)/V3(19» 17) /VB0{19»17 )» P3( 19/17) / PBO (19/17 
i)/W6(19/17)/PBAa9/17)/WR0(900)/ AN{15)/BN(15)/CN(15)/CS( 191/KAP (17 
S)/Z{17)/CSA(19)/DCA(17)/DR( 19 »/CORIOL (19 )/ GAM8< 19) / FM( 19/17)/FT(19 
$/ 17)/ 08(19/ 17)/ DcNS(17)/CNA(19)/CNBC19)/CUB(19)/TN(19)/DCDS(19)/ 
$RAYF (17)/CH1(19/17)/RR (13/15)/ AMU8)/ 6t1(16)/C)Uia)/GHEP(19)/XeA(l9 
$/17)/ WBD( 19/17)/GBV(17)/DUH(19) 

REAL PHB1(19}/XP,1C19)/0EL(19)/$TAB(19)/ALA(19)/ALB (19)/TNA< 19) / $U 
$a9)/SV(19)/GAtU19)/CDRa9)/XMAl(19)/NGC (19 ) / A1 ( 17 )/C 1 ( 17)/ ANG( 19 ) 
S/TNB (19)/WRA( AOO)/AR (15)/ BR ( 15 ) / CR (15 ) / T6 ( 1 9, 17 ) 

REAL KZ/TZ0(17) 

INTEGER IPHAS(19)/S/SYM 

COMMON UB/PB/UI/UIO/ VI/ VIO/ Pl/PiO/ Wi/ UBO/ V80/ PBO/WB/QB 
COMMON /BUMLEG/ CURTIS { 13/ 19 ) / T6 ( 16 )/ T ( 1 B/ 1 9 > / PRSG ( 1 8 ) / ? RS ( 18/ 19 ) / 
$QDG{18)/Q0(18,19)/Q0ZS6(1B)/Q0ZS ( 18/ 19 ) /CTS ( 18/ 19 ) / CG ( 1 8 / 19 ) / ZN (1 8 
S)/ AZN(18) /FDY(18)/TQ(18/16)/CZT(19/17),SZT(19/17)/ CDT( 19/ 17 )/ RED( 6 
$/16)/0RS(16) 

DIMENSION QA(17) 

EQUIVALENCE (XBA/PBA)> (TB/PBA) 

S » 1 

SYH -0 FOR global DOMAIN/ SYH "1 FOR ANTISYMMETRIC HEMISPHERE 
SYM » 0. 

DATA JM/KN/ IF LG719/ 17/ 0/» RAD/ D2/ BVF/SH/& . 37c6/ 5000 ./ A. E-A/7 . 0E3/ 
CELT • IBOD.AZ. 

DT - DtLT/2, 


IFD « 46 
ICT “ IFO-1 
KZ - 0. 

SETUP CONSTANTS AND INITIALIZE 
STRDAY • 80. 


READ 

(2/196) 

CURTIS 


READ 

(2/200) 

T6/T 


RFAD 

(2/160) 

({PRSd/ J)/I»l/4)/ J 

=1/19) 

READ 

(2/195) 

C3G/Q0 


READ 

(2/195) 

QOZSG/ODZS 


READ 

(2//.95) 

CTS/CG 


READ 

(2/205) 

RED 


WRITE 

(6/185) 

{ (REDd/ J)/ I«I/6)/ 

J-1/16 


DATA INIT/IEND/ IPRINT/1/ 1/1/ 
PI « 2,^AS1N(1. ) 

JML " JK-1 
KNL « KN-1 


M « 
N ■ 
DY • 
Cl > 
CZ2 
PY2 


JK-2 

KN-2 

PI*RAD/JML 

(O./l.) 

■ 0Z**2 
» 0Y**2 


FREQ » 7.292E-5 
RHZ - SH/287,1 
DTOOr » DT/{9.*DY) 

EPL « exP(DZ/(4.*SH) ) 

EMI » EXP (-DZ/< 4.*SH) ) 

DO 10 J-l/JM 

SB - PI*(JH-2.*J) /(2,»JML) 

SA » PI*< JM+1-2.*J )/(2.*JML) 

CS(J) » C0S(S8) 

CSA(J) « CDS(SA) 

CORN) - 2.*FREQ+SIN(SB) 

CORIJL(J) - 2,*FP60*SIN(SA) 

GAKJ) - l./( l.+CDR(J)**2*0T**2) 

CN3(J) » l./( DY4CS( J)) 

PH8*(J) - 150.f (51N(3.«(SA-PI/6) ) )**249.8 
IF (SA.GT.-PI/6.) PHBl(J) » 0. 

ALA(J) « DT/DY*COR{J) 

AL3(J) ■ COR( J )*DT*GAM(J) 

XMl(J) « S/(4A0*CS(J ) ) 

TUlU) “ GAM( J )*CI*XM1( J)+DT/2, 

TVl(J) « TUK J)*COP(J)*DT- 

GAMBU) • l./U.+C0RI0L(J)»42*0T**2) 


PAGE KS 
QUAinrY. 


A 1 
A 2 
A 3 
A 4 
A 5 
A 6 
A 7 
A 8 
A 9 
A 10 
A 11 
A 12 
A 13 
A 14 
A 15 
A 16 
A 17 
A IS 
A 19 
A 20 
A 21 
A 22 
A 23 
A 24 
A 25 
A 26 
A 27 
A 28 
A 29 
A 30 
A 31 
A 32 
A 33 
A 34 
A 35 
A 36 
A 37 
A 36 
A 39 
A 40 
A 41 
A 42 
A 43 
A 44 
A 45 
A 46 
A 47 
A 49 
A 49 
A 50 
A 51 
A 52 
A 53 
A 54 
A 55 
A 56 
A 57 
A 56 
A 59 
A 60 
A ei 
A 62 
A 63 
A 64 
A 65 
A 66 
A 67 
A 68 
A 69 
A 70 
A 71 
A 72 
A 73 
A 74 
A 75 
A 76 



CNA(J) • BVF*DZ240T/(DY*CS(Jn 
CUB (J) - GAMB(JJ*CSACJ) 

SV(J) » GAH(J)*DT/DY 
SU(J) « COR(J )*DT*SV(J) 

K6CU» ■ GAK( J)*CS(J) 

TNA(J) - (CSA(J)/C5(JJ-CS{J)/CSA(J))/0Y 
10 NUl(J) « J)*GAH( JJ^CSfJ) 

DO X3 J-2>JM 

DEL(J) - BVF*DT*DZ2;CSA(J) 

STA6(J) * DEL(J)*DT 
DCOS(J) “ l./(RAD*CSA(Jl) 

XHAKJ} » S/( RAD*CSA<J) ) 

6HKJJ • .S^CIAXHAH J)*DT 

TNB(J-I) » {CS(J-1>/CSA (J J-CSA(J)/CS{ J-1) > /OY 
15 TK(4» « l./C4.*DY)*(CSA{J-l)/CSA(J)-CSA(J)/CSAU-in 

DO 20 K»1>KN 

Z(K) - (K-1JA=DZ 

D£NS(KJ » EXP(Z(K)/<2.*SH) ) 

6BV(K) B 1. 

IF (K.GT.KML) GO TO 20 
KAP(K) « £RED(5,K)/e6-t00.)*DT 
20 CONTINUE 
IRAD ■ 0 
IRCT - 1 

WRITE {6»I<35) {KAP(K),K-1^KN} 

DO 25 J«l,JhL 

DZN « BVF*'DZ2/<CS (J )*DY2)*DTt'*2 
AHU) “ DZN*GAMB (J')’tCSA(J) ' 

CH(J) • 0ZN*GAHB(J+1)*CSA( J+1) 

25 BH(J) - -AHtJJ-CKJJ 
BH(1) - BH(1)+A«a> 

BH(JHL) - BH( JML)+CM< JML) 

AHdJ « CH(JML) “ <0,,0. ) 

ITIHc B 0 
IKAT • 0 
HT • 1. 

TIME - 0, 

IF (INIT«NE.O) 6Q TO 90 
QRS<I» « -XOI. 

C READ INITIAL ZONAL FLOW 
DD 30 

Y - 2.«PI*STR0AY/360. 

RZ - ?RS(1/J3 
R1 - PRS«2,J) 

SI « PRS(3#J3 
R2 « PRS(4>J} 

UBOU/1) « RZ/2,+Rl*C0S{n+Sl+SIN(Y>+K2+C0S(2.+Y) 

30 CONTINUE 
DO 4J 

DO 35 K-1>KN 

UBO(JiK) - UB0(J#1) 

B UBOIJ^K) 

35 continue 
AO continue 

C CCHPUTE INITIAL GEOPOTENTIAL 
DO 55 K-1»KN ' 

DO A5 

U80{J»Kl - U3D{J# K)/DENS(K) 

45 . CONTINUE 

PBA(X,K) « 0. 

00 50 J«2»JHL 

PBA(JjK) « PBA(J-X>K)+DY*tCORIOL( J)*UB( J>K)-UB (J»K)*(U8t 
$ ♦TN{J)+UB( J+X»K) ♦TNI 4+1 n ) / dens (K) 

50 CONTINUE 

55 CONTINUE 

DO 60 K-1»KN 

PBA(JM,K> B P8A(J«L»KJ 
60 CONTINUE 

DU 75 K-1>KN 
SUMi - 0, 

TEMP • 0. 

DO 65 JbX^JNL 

TEMP • T£f.°+CS(J) 

SUMi ■ SUMi+PBAU^KI+CS { J) 

65 CONTINUE 

00 70 4«1/4KL 


A 77 
A 78 
A 79 
A eo 
A 61 
A 62 
A 83 
A 8A 
A C5 
A 66 
A 87 
A BB' 
A 69 
A 90 
A 91 
A <52 
A 93 
A 94 
A 95 
A 95 
A 97 
A 98 
A 99 
A 100 
A 101 
A 102 
A 103 
A 104 
A 105 
A 106 
A 107 
A 108 
A 109 
A 110 
A 111 
A 112 
A 113 
A 114 
A 115 
A 116 
A 117 
A 116 
A 119 
A 120 
A 121 
A 122 
A 123 
A 1Z4 
A 125 
A 126 
A 127 
A 126 
A 129 
A 130 
A 131 
A 132 
A 133 
A 134 
A 135 
A 136 
A 137 
A 138 
A 139 
A 140 
A 141 
A 142 
A 143 
A 144 
A 145 
A 146 
A 147 
A 148 
A 149 
A 150 
A 151 
A 152 



PBA(J,K) - P3A(J»K)~SUH1/TEHP A 1»;3 

70 CONTINUE A 154 

P8A(JH*K) • PBAUHL#K> A 155 

75 CONTINUE A 156 

DO 85 J-1,JM A 157 

DO 80 K-1»KN A 158 

P8<J»KJ « PBACJ^K) A 159 

PBOtJ>K) « PBA(J^K) A 160 

80 CONTINUE A 151 

85 CONTINUE A 162 

GO TO 95 A 163 

C A 164 

90 CONTINUE A 165 

C FOR CONTINUATION RUNS READ INPUT DATA A 166 

READ (I) TIMEf Pij.PB^Pl0,Ul^UB/U10,Vl^ V1DW1>PBO/UBO,VB^VBD>WB»WBO, A 167 

SGBV/TZOjZN, AZN, FDY»TQ»ORS>SZT A 168 

95 CONTINUE A 169 

C GLOBAL HEAN STABILITY PROFILE AND INITIAL HEATING A 170 

OAU) = 101. A 171 

CALL HEAT ( TIME# GBVr PB> TB#QB»RHZ, OENS»DZ» J R> KN> TZQ» 8VF/ SH* PI,OA, A 172 

SSTRDAY) A 173 

DO 100 K-l/N A 174 

ANIK} =• GBV(K) A 175 

CN(K) « G8V(K+1) A 176 

BN(K) - -66V(K)*EPL**2-GBVCK+1J*EKI**2 A 177 

100 CONTINUE A 178 

DO 105 K«1,N A 179 

ARdO ■ AN{K) A leO 

CR(KJ « CN(K) A 161 

105 6R(K) • BNIK) A 1E2 

BR(N) » SR{N)+EHI/EPL*GBV{KNU A 163 

ARID ■ CRIN) » 0. A 184 

DO 110 J»2,JML A 185 

Al(J-l) - STASU)’»‘(NGC(J-i)/DY2>-XMi{J-l)**2+NGC(J~l)/4.) A 1B6 

ClfJ-1) - STAa(J)+(NGC(JJ/DY2-XKl{J)**2*NGC(J>/4.) A 187 

Bl(J-l) - -5T AB( J) + ((XH1( J-Dt + E^NGCC J-1 ) + XHl( J)+*2*NGC (0 ) )/4.+ A 188 

i (N6C(J-1)+N6C(J ) )/Dr2+Cl*(KGC(J-l)*XNi(J-lJ+ALAU-l)-NGC(J )*Xni A 189 

i {J1*ALA(J))) A 190 

110 CONTINUE A 191 

C CDNPUTE EODY NONNEtfTONION HEATING A 192 

115 DO 120 A 193 

DO 120 K*1>KnL a 199 

120 QS(J»K) « I0.#0.) A 195 

DAY « TIHE/(24.*60,*60.) A 196 

DO 125 J«1»JH A 197 

y - 2,*PI*{STRDAY+DAY+DELTUe64o0.»/360. A 193 

RZ • PRSCl^J) A 199 

R1 » PRS(2^JJ A 

Si • PRS(3»U) A 201 

R2 • PRS(4»J) A 202 

PR5(5#J) » RZ/2,+Rl*C0S (Y»+Sl*SIN<Y)+R2*CaS(2«*Y) A 203 

125 CONTINUE A 2U4 

PRS{6^i) ■ 0. A 205 

DO 130 J»2#JML A '206 

PRS(6»J) - PRS(6/J-l)+DY*C0RI0L(J)*PRS{5#J)-PRS(5jJ)*(PRS{5»J-l) A 207 

5 ♦TKfJ)+PRS<5> J4-l)ifTN(J+ll) A 20B 

130 CONTINUE A 209 

PRS(6,JH) - PRS(6#JML> A 210 

SOM » 0. A 211 

TEMP « 0, A 212 

DO 135 J»1#JKL A 2J3 

TEMP ■ TEHP+CS(J» A 214 

SUM • SUM+PRS(6/J)+CSU) A 215 

135 continue a 216 

DO 190 J-l»Jf1L A 217 

PRSC6>J) « PRSC6> J)-SUM/TEHP A 218 

140 CONTINUE A 219 

PRS<6»JM) « PRSt6>JML) A 220 

00 145 A 221 

PRS(7»J) « (PRS(6^J)+2.*PB(J,1)+PB0(J»1) J/4, A 222 

145 CONTINUE A 223 

DC 15CJ K-1>KN A 224 

RAYFdO - DT*(O.5/(96O.*36U!j.)+{l.+TANHnZ(K)-55w00.)/10u00.))/ A 22 = 

$ (96,*3600.))*(l.-EXP{-.4£-05*TIri£)) A 226 

15U OCA<KJ * K4P(K)»(1,-EXP(-(TIHE+DELT)*,4£-5)J A 227 

DKY - I.£-b*(l,-cXP (-,4£-5*TIM£) ) A 226 



DELAY * 0. 

DELTIH - TIME-DELAY 
GTIHE - A.32E+05 
IF (DELTIM.LT.O.) 6Q TO 16Q 
DO 155 

PLKJ) • (PHBHJ )*(1.-EXP(-(DELTIH+DELT)/GTIHE))+2.*P1(J#1) + P10 
S 

155 CONTINUE 
160 CONTINUE 

DO 165 K-l/KNL 

P1A(1»K) « tO.fO.J 
165 PlAUKfK) » (0,^0.) 

CALL EDDY ( Al, 61, Cl» ALB> B VFj» ANG# Cl, CS^ C SA> COR» DCA» DELjDSNS» OKY> 
$DTODY/OT,OZ,DY,EMI,£PL,GAr.,GHl,ICT, iFLG, I FD, J H, KNL, H# S, NGC, NUi, 
$RAYF>PLl,P»0,SJ»SV,TUi,TVl,TNA,TNB,XMl, XMAl,PB,PlA,Pl,PlD»R,UB#Ui» 
SUIO, V1,V10,W1,*<R0 ,CDUH,AN,BNj.CN,GBV,SYM,QS,KZ ) 

CALL FLUX ( JM,KN,DV,DZ,CSA,CS,D£NS, cPL, cKI, Ul, VI, Wl, PI, Pi A, FH, F T, 
$KZ,UBO,COUH,CNB) 

IFLG « IFLG-2 
ICT - ICT-1 

CALL ASTREAM ( AM, B M, CM, AR, BR, CR, IFLG, IFD, ICT#DT,OZ, DY, CHI ,QB,WRA, 
$RR,XBA,UB,UBO, VB, VBO,WB,WEO,P B,P BO, PBA, CS, C3 A, FH, FT, E PL, EHI,RAYF, 
tOCA,DKY,TN,DR,DENS,CNB,6KtP,CORIOL,CNA,CUB,BVF, JH, JML,KN,N,GBV) 
IHAT » lHAT+1 
ITIHE » ITIME+1 
TIME « TIME+OELT 
IF (ICT.EO.IFO) ICT » 0 


170 


C 

175 


C 

c 

leo 

165 

1<?0 

195 

200 

205 


IRAD « IRAD+1 

IF (XRADtLT.IRCT) GO TO 170 
OAU) - 0. 

CALL HEAT ( TIME, G6 V,PB, TB, QB, RHZ, DENS, DZ, JK, KN, TZO, BVF, SH, PI, QA, 
$STRDAY) 

IRAD « 0 
CONTINUE 

IF (IMAT.EQ 4 IPRINT) GO TO 175 
GO TO 115 

CONTINUE 

CALL AOUT (TI‘'E,DT,KN, JH,DUK,D£NS,IPHAS,Z, JML,0Z,RHZ,C0R,DY,CI,XH1 
$,IENO, 1TIME,MT, IMAT,P1, PB, Ul, UB, Vi, V3, Wl, Ufo, FH, FT,S,QS,CS,Q3,CZT, 
$SZT,COT,Ta) 

REWIND 1 

WRITE (H TIM£,Pl,PB,P10,Ui,UB,UlU,Vl,V10f *1, PBO,l)BO,V3,VBO,Wfi,WBO 
S,GBV,TZO, ZN,AZN,FDY,TO, QRS,SZT 
END FILE i 
END FILE 3 

IF (ITiMt.LT.IEND) GO TO 115 
STOP 


FORMAT {AF15.A) 

FORMAT (5X,6F1A.7) 

FORMAT (oF10« A, /,5F1C.A) 
FORMAT (9E10.3) 

FORMAT {9F10.3) 

FORMAT (6F10.0) 

END 


A 229 
A 230 
A 231 
A 232 
A 233 
A 23A 
A 235 
A 236 
A 237 
A 238 
A 239 
A 2A0 
A 2A1- 
A 2A2 
A 2A3 
A 2AA 
A 2A5 
A 2A6 
A 2A7 
A 2A6 
A 299 
A 210 
A 251 
A 252 
A 253 
A 259 
A 255 
A 256 
A 257 
A 258 
A 259 
A 26Q 
A 261 
A 262 
A 263 
A 264 
A 265 
A 266 
A 267 
A 2t« 
A 269 
A 270 
A 271 
A 2 72 
A 273 
A 274 
A 275 
A 276 
A 277 
A 278 
A 279 
A ZFO 
A 2tl 
A 2B2 
A 263 
A 2t4 
A 285 



c 


c 


c 

c 

c 


10 


c 

15 


C 

20 

25 


30 

C 


35 

^0 

^5 


50 

55 

6C 

C 


SUBROUTINE HE AT {TIHE»6SV» PB»TB> QB^ RHZ» 0 ENS/ 02, JM#Kh>TZO> BVF,SH, PI> 
SQA,STRDAY) 

COMKON /BUMLEG/ R (13, 19 » , TGI 18) >TC Id, 19), PRSG {18 ), PRS (13, 19 ) , QOG 
S(16),00(18,19),QQ2SG(lS),ClDZSa6,19),CS{18,19), C6(16,19),ZN(le), 
5AZN(18), FDYtl8),T0(ie,16),CZTC19,17»,S2m9, XT) , COT (19, 17 ), RED ( 6, 
$16),QRS(16> 


DIMENSION GBV(KN), PB(JH»KN), TBtJK,KN)> QB(JM,KN), DENS(KN), TZO 
$(KN), QAIKN) 

JML » JM-1 
KNL « KH“1. 

LBOT - A 
tTOP » 19 

COMPUTE SOLAR SEOHETRY FACTORS 
DAY - TIME/(2A.*60.*60.) 
lOAY “ DAY 
RES « DAY-IDAY 
IF (RES*GT,0,05) go to 60 
PERH = (101,21972+180, >*360. /360.+80* 

VD « STRDAY+DAY-PERH 
V » V0*2»*PI/360. 

EC - 0.016722 

RHO « (l.-EC*£C)/ (l.+£C*CDS(V) ) 

DELTA-SOLAR DECLINATION 




DELTA » O,A093198*SIN(2,*PI*(STR0AY+DAY-8O.)/360.) 

DO 30 L-1,18 
LL » 19-L 
RL.- LL-9 
PHD « RL+10.-5, 

PHI-TERRESTRIAL LATITUDE 
PHI « PHD+PI/ISO. 

TSTAR-TIHE OF SUNSET (OR NEGATIVE TIME OF SUNRISE) 
ZEN-AVERAGE VALUE OF COS(SZA) BETWEEN -TSTAR AND TSTAR 
SUN - TAN(DELTA)*TAN(PHI) 

ISUN - SUN 
IF (iSUN) 10,15,20 
TSTAR - 0, 


ZEN - COS(DELTA-PHI) 
GO TO 25 


TSTAR - (12./PI )*ACCS (-SUN) 

ZcN - S INC DELTA) *SIN (PHI ) + (12. / (PI*T3TAR) )*CQS( DELTA) *005 (PHI )* 
$ SIN(PI*TSTAR/12.) 

GO TO 25 

TSTAR - 12. 

Z£N - SIN(DELTA) *S1N(PHI) 

CONTINUE 

AZEN « 33,/S0RT(1224.*Z£N*ZEN+l.) 

FDAY - TSTAR/12. 

ZN(L) - ZEN 

AZN(U » AZEN 
FDY(L) - FOAY 
CONTINUE 

COMPUTE GLOBAL RADIATIVE EQUILIBRIUM TEMPERATURE 
IF (TIHE.GT.O.) GO TO A5 
DO AO K»1,KNL 
TZO(K) - 0. 

DO 35 J-i,JKL 

THETA - 175-(J-1)*10 
THETA - THETA*P1/180. 

T2000 « TZ0(K)+TU,K+3)*SIN(THETA)*PI/38. 

CONTINUE ^ 

CONTINUE 

CONTINUE \ 

DO 55 J=1,JM ^ 

DO 50 K-1,KNL 

SZT(J,K) « S0L(TZQ,JpK+3,RH0) 

CONTINUE 
SZT(J,!LN) • 0. 

CONTINUE 

CALL RAOEQU (TZO, TB, CB, RHO ) 

CONTINUE 

CC,’-?UTE ZO'-AL MEAN TcMPERaTURc (DEVIATION FROM ‘GLOBa£. AVERAGE) 

DO 65 J-1,JKL I • 

DO 65 K-1,KNL ' 

TB(J,K) - RHZ* (P3{J,K+1)*D£NS(K+1)-PB( J,K)'*OENS(K) )/0Z 


B 1 
B 2 
B 3 
B A 
B 5 
B 6 
B 7 
6 8 
B 9 
B ID 
B 11 
.B 12 
B 13 
B 14 
B 15 
B 16 
B 17 
B 18 
B 19 
B 20 
B 21 

5 22 

B 23 
B 24 
B 25 
B 26 
B 27 
8 28 
8 29 

B 30 
B 31 
B 32 
B 33 
B 34 
B 35 
e 36 
B 37 
B 36 
B 39 
B 40 
B 41 
B 42 
B 43 
B 44 
8 45 

8 46 

P 47 
B 48 
8 49 

B 50 
B 51 
B 52 
B 53 
B 54 

6 55 

8 56 

B 57 
B 58 
B 59 
B 60 
B 61 
B 62 
B 63 
B 64 
B 65 
B 66 
B 67 
B 6E 
B 69 
B 70 
B 71 
B 72 
e 73 
8 74 

B 75 
B 76 


65 




IF ITIHE.CT.O.) CO TO 75 

6 

77 

c 

COHPUTE SUTtC STAEUITY FXOH TZO 

B 

78 


H • KNL-1 

B 

79 


DO 70 K«a»N 

B 

60 

70 

66V{K) « SVF*RHZ/t2./7*nZ0(Kl/SH+(TZ0<K+l)-TZ0(K-l))/l2,*DZn 

B 

81 


6BVa) > 66V(2) 

6 

82 


G6V<KKL) > GBVCN} 

6 

63 

75 

continue 

B 

BA 


IF <QA(1).LT*100. ) GO TO 80 

B 

85 


WRITE (6>105I (GBV(I)#I-1»KNL) 

B 

B6 

60 

CONTINUE 

6 

67 

C 

CONFUTE HEATING RATE 

6 

86 


CALL OOOT (TZ0»TB/Q6f RKO) 

6 

69 


GTIMc - 1.73E+06 

> B 

90 

• 

DO 90 K«1»KNL 

B 

91 


DO 65 L-1»JNL 

6 

92 


QB(L/K) • OB(L»K>/CB6AOO.*RHZ*DENS<K>>*I1.-EXP{-TIME/GTIK£))* 

B 

93 


» DZ 

B 

9A 

85 

CONTINUE 

B 

95 

90 

CONTINUE 

B 

96 


OAtl) ■ 0. 

B 

97 


OA(IS) • 0« 

B 

96 


aA(16) * 0. 

B 

99 


DO 100 H-L80T»LT0P 

B 

100 


I » N-3 

B 

101 


QACI) - 0* 

B 

102 


DO 95 L«1#JML 

B 

103 


THETA - 175“IL-1>*10 

B 

104 


THETA ■ THETA+PI/ieO. 

B 

105 


QMI> ■ QA(I>+QB(Lf I)*SIN{THETA)*PI/36, 

B 

106 

95 

CONTINUE 

B 

107 

100 

CONTINUE 

B 

108 


RETURN 

B 

109 

C 


B 

liO 

C 


B 

111 

105 

FORMAT (5X»eE15,3) 

B 

112 


END 

B 

113- 



FUKCTIQN DELT(TP,L#IH>RHO» 

.COMMON /BUHLEG/ R ( 13» 19 ) ,TG U8 J »T( 16> 1<5) , PRSG (16 )# PRS ( 13» 19 ) »QOG 
sa8)>QQll3,19)»Q0ZSG{18),Q0ZStX8,19)iCS Cl 8# 1 9 ) , CG ( 16» 19 ) , ZN ( 16 ) > 
SAZK(16WFDY(lSWTQ(18»16)»CZm9#17J>SZT(19>17)>CDT(19,17)»RED(6> 
$16)>0RS(16) 

DIMENSION TP(17)» TDC19) 

I - IH~3 
COZONE » 0, 

CCD2 - 0. 

S03 - 0, 

IF (I.EQ.16) GO TO 20 
KNL " 16 
DO 10 4»1,KNL 
T0CJ+3J *= TPCJ> 

10 CONTINUE 

C FOR Z .LT 18.5 KM U, S. STANDARD ATMOSPHERE ASSUMED 
DO 15 J=l#3 

TD(JI T(L,J) 

15 CONTINUE 

SOS - SZT(L^I3 

C DICKINSON INFRARED COOLING 
TR - RED(3,I) 

OR • ORSd) 

A « R£D(5»1) 

CCD2 * -QR-A*(TD(1M)-TR> 

20 CONTINUE 

CDT(L<n «• CC02 

CZT(Ld) * COZQNE 

DELT « S03+CCQ2+CDZDNE 

RETURN 

ENTRY SOL 

I • IH-3 

OZ - 5000. 

SH - 7000. 

Z * a-iJ*DZ 

AIRDEN « 1.539l£-04+EXP(-Z/$H) 

ZEN • ZN(L) 

AZEN » AZN<L) 

FDAY » FDY(L) 

TAU » 10. 

CLOUD » 0.446 
R6 ■ 0,5 
ICLQUD ■ 2 

C HEATING BY THE ABSORPTION OF SOLAR RADIATION BY OZONE 

C FROM LACIS AND HANSEN, J ATMOS SCI 31 116-133 

S03 « 0, 

IF (FDAY.LT.O.OOOIJ go to 25 
C CLEAR SKY 

RAB « 0.219/(1, +0.6164ZENJ 
ROB • 0.144 

RBI - RA3+(l,-RA9>»(l.-RDfl)*RG/(l.-RDB*RG) 

RB - RBI 

C CLOUDY SKY 

RAB - O.I3*TAU/ (l.+0,13*TAU) 

POB « RAB 

RBI » RAB+(i.-RA3)*(l.-RDB)*R6/(l.-RDB*RG) 

UT • AZEN*QOZS(L,ICLOUD»+1.9’M002S(L,ICLQUO>-QOZS<L,IM)) 

A1 - QZUV(UT} 

UT - AZ£N*0aZS3<L)+1.9*(00ZSG(U-QG2S(L,IK)> 

A2 - OZUV(UT) 

U - AZEH+QQZS(L,IH) 

A3 - OZUVCUl j 

S03 « 76.374*Q0(L,IH)n.t+05*(l./AIRDENl*ZEN*FDAY*(AZEN*A3'Kl,- 
$CL0UD)*RB+1.9+A2+CLaUD*RBl*1.9*Al) 

S03 ” SQ3/(22G.*2.67*PH3#RH0) 

25 CONTINUE 

DELT « SQ3 

RETURN 

END 


C . 1 

C 2 
C 3 
C 4 
C 5 
C 6 
C 7 
C 8 
C 9 
C 10 
C 11 
C 12 
C 13 
C 14 
C 15 
C 16 
C 17 
C 18 
C 19 
C 20 
C 21 
C 22 
C 23 
C 24 
C 25 
C 26 
C 27 
C 26 
C 29 
C 30 
C 31 
C 32 
C 33 
C 34 
C 35 
C 36 
C 37 
C 38 
C 39 
C 40 
C 41 
C 42 
C 43 
C 44 
C 45 
C 46 
C 47 
C 48 
C 49 
C iO 
C 51 
C *2 
C 53 
C 54 
C 55 
C 56 
C 57 
C 53 
C 59 
C 60 
C 61 
C 62 
C 63 
C 64 
C £5 
C 66 
C 67 
C 68 
C 69- 



FUNCTION 02UV(UJ 

SUBROUTINE TO CALCULATE SPLAR ENERGY ABSORBED BY OZONE 

FI - 1.0+(O.OA2tUJ + (3,.23E-4*<U**2.0n 

FE - {O.0E12/Fl)*{1.0-((U/Fl)l‘(O.0'*2+(6.46E-4*U)n ) 

FI - 1.0+(13b.6*U) 

F2 « F2+{1.082/(Fl**0.eo5))*a,0-(a38.6*0.805*U)/Fl)) 
FI » 1.0+( aos.b^^ui^^s^o) 

OZUV - F2+({0.0656/Fl) + {1,0-<3.0+(FI“1,U)/Fi)n 

RETURN 

END 



o o 


SUBROUTISH RADEQUtTZDiTBj QB jRHO) E 1 

DIKENSIQN TZ0(17)/ TP(17)j TBtl9,17)» QB (19/ 17)> ' A3{ 2> E 2 

COMMON /SUHLEG/ R ( 13/ 19) > TGUe ), T (1 5/ 19 ), PRSG C18 )/ PRS 1 18j 19 >, QOG E 3 

J(18 »*QQa8/19)/0aZSG(18)»Q0ZSU8»l9)/CTS (18/ 19), CG (18/ 19 )»ZN (16 >i E 

$A2Na6)/F0Y{lB)/T0{18/16)/CZT(19/17)/SZT<19/i7),CDT(19/17)/RED(6/ E 5 

$16)/QRS(16) E 6 

DATA LBQT/LTDP/ ITOT/DT/EPS/4/19/20/0. 1,0,1/ E 7 

DATA JML/KNL/16,16/ E 8 

PI ■ AC0S(-1«) E 9 

IF (QRS(l).GT.-lOO. > GO TO 20 E 10 

DO 15 K-L80T/LTQP E ll" 

I « K-3 E 12 

QRS(I) -0, E 13 

DO 10 J-1,JKL E lA 

THETA » 175~(J-1)*10 E 15 

THETA •= THETA*PI/180. E 16 

W » SIN(TH£TA)«PI/36,’ E 17 

QRS(I) » QRSCI)+S0L(TZ0/J/K/1.)*W E 16 

10 CONTINUE £ 19 

15 CONTINUE E 20 

20 CONTINUE E 21 

DO AO IT-l/ITOT E 22 

SUM * 0. E 23 

LHP - LTOP-1 E 2A' 

DO 35 H-LBDT/LMP E 25 

I « H-3 E 26 

TZO(I) - TZ0(I)+2,*DT E 27 

DO 30 J«X/2 £ 28 

TZO(.I) - TZO(I)-DT E 29 

A3(J) - 0. E 30 

DO 25 L«1/JML E 31 

THETA - 175-(L-1)*10 E 32 

THETA - THfcTA«PI/100, £ 33 

U - SIN(THETA)*PI/36. E 34 

A3(J) ■ A3( J)+DELT{TZO/L/M/RHO)+W E 35 

25 CDUTINUE E 36 

30 CONTINUE E 37 

DQ " (A3(l)-A3(2) )/DT E 36 

OIFF e A3(2)/DQ E 39 

TK - TZOm-DiFF E 40 

DTP - OIFF £ 41 

IF (0TP.GT.20.) DTP * 20. E 42 

IF (0TP.LT.-20. ) OTP « -20. E 43 

TN ■ TZO(n-DTP E 44 

OIFF » AflS(DIFF) E 45 

IF (D1FF.GT.3UM) SUM - DIFF E 46 

TZO(I) - TN E 47 

35 CONTINUE E 48 

IF (SUM.LT.EPS) GO TO 45 E 49 

40 CONTINUE E £0 

PRIM 65/ ITOT E 51 

STOP 1 E 52 

45 print 70/ IT . E 53 

WRITE (6/75) TZO E S'* 

RETURN E 55 

ENTRY QDQT E 56 

C COMPUTE RADIATIVE HEATING IN KELVIN PER DAY E 57 

00 60 L-l/JML E 58 

DQ 50 J-l/KNL E 59 

TP(J) ■ TB(L/J»+TZ0U) E 60 

50 CONTINUE E 61 

0B(L/1) * 0. . E 62 

DO 55 M-LBDT/LTQP E 63 

I - H-3 E t** 

QB(L/I) » DELKTP/L/M/RHO) E 65 

55 CONTINUE E 66 

60 CONTINUE E ^>7 

RETURN E 68 

E 69 

E 70 

65 FORMAT (5X/+TEMPERATURE PROFILE FAILED TO CONVERGE AFTER*/I3/ E 71 

i* IIERATIONSi- ) E 72 

70 format ( 5A/*TtMP£RATURE CONVERGED AFTER */I2/* ITERATIONS*) E /3 

75 format (iX/18F7,2) E 74 

end E 75- 



o o o 


SUBROUTINe ASTREAM(iM#BM»CK»AN»BN»CN* IFLG,IFO»lCT,DT»OZ»OY»CHI»<IB, 
SWRO»RR»XBA»UB»U3Q^VB/VaO/K3/WBO»PBjPBCI»PBA*CS>CSA>FM»FT*EPL^EHl>- 
tRAyF/DCA>OKY>TN,DR,DEHS»CNB»G«EP>CURIOL>CHA>CU3#BVF# JM»JML»KN>N» 
$6BV) 

C 

• C 

OIHEKSION AM(JHL)# CfUJMD^ PBA(JK>KN)» CQRIOL(JH)> QB(JH 

$>KN), OCA(KN>» UB(JH^KN), UBO(JM^KN)> VBtJrtiKM)# VBD(JM>KN)> PB(JM 
SfKNJ/ PBO(JM/KN>» TN(JH)» DR<Jri)> CHI{Jh>KN)# WRO(AOO)> 

i RR(JML»15), GM£P{JH)/ XBA(JH>KN)» OENS(KN)/ CS(JHJ» C5ACJH), FH 
S(JH»KW)# FT(JH»KN)» CNA(JH># CNB(JH)> CUfi(JH>i RAYF(KN)# MBOUHjKN 
SJ/ A«a5)# BN{15)» CN(15)» GBV(KM) 

COMMON /BUHLEG/ CURTIS(13,19)#T6a8J»T(i3#19)#PRSG(i8),PRS(l8»19) 

M » JMl-1 
KNL « N+1 
KH - N+2 
C 

C CHOOSE LEAPFROG OR FORWARD DIFFERENCE 
C 


DO 10 J“1>4M 
DO 10 K-2,KN 

10’ U8{J>K) - UB tJ»K) /DENStK) 

ICT • ICT+l 

IF (ICT.LT.IFDJ GO TO 15 
T1 » 1. 

T2 « 0. 

T3 - 1.- 
TA « 0, 

T5 « 2. 

GO TO 20 
C 

15 TI » T2 « .5 
T3 - 2. 

TA - 1, 

T5 « 4, 

20 CONTINUE 
C 

C UB SMOOTHING 






F 1 
F 2 
F 3 
F A 
F 5 
F 6 
F 7 
F 8 
F 9 
F 10 
F 11 ■ 
F 12 
F 13 
F 14 
F 15 
F 16 
F 17 
F 16 
F 19 
F 20 
F 21 
F 22 
F 23 
F 2A 
F 25 
F 26 
F 27 
F 28 
F 29 
F 30 
F 31 
F 32 
F 33 
F 3A 
F 35 
F 36 
F 37 
F 38 


C F 39 

DO 30 K“2jKNL F AO 

DO 25 J«3jH F A1 

25 « FH(J,K)-DKY*(U30lJ-2,K)/CSAU-2)-A.*UB0{J-l,n/CSACJ-l F A2 

$ ) + 6.*UB0C J»K) /CSA(J )-A.*UBO(J + l» K)/CSA{ J+lJ+U80(J+2»K)/CSAU+2) ) F A3 

$ /CSA{J)«+2 F AA 

FM(3»R) « FMO^Kj-DKY+tUBOtZ^Kj/CSAta) )/CSA (3)**2 F A5 

FK{K#K) - FMMf lO-*DKY+tUBO(4rtL,K)/C3A ( JHU) /CSA(M)**2 F A6 

FM(2*K) • FM(2»RJ-DKY+( 2,*UBa(2#K)/CSA(2)-3t*UB0(3»K)/CSA{3)+UBD F A7 

S (4»KJ /CSA(A)) /C5AC2 >**2 F 49 

FM(JKL#K> o FK( JHL»K)~DKY*(UBa(JKL-2»K)/CSA{JHL-2)-3,*UB0(JHL-i/ F A9 

$ K)/CSA(JKL-1)+2.*UB0(JHL^K)/C5AIJHU )/CSAUHL)**2 F 50 

30 CONTINUE ' F £1 

F 52 

THICKNESS TENDENCY F 53 

F £A 

DO 40 F 55 

DO 35 K«1»KNL F 56 

35 PBA(J»K) « F8( J,K+1)*EPL-PB(J»K)*EMI F 57 

40 PBAU^KN) - 0. F 58 

CO 65 K-i»KNL F 59 

DO 50 J«1^JML F 60 

CHi(J>K) ■ Tl’f'PBAC J>K.)+T2*(PBO(JjK+i)*tPl-PBO(J»Kl*EMI)+DT*(QB F 61 

i C J>K)+FT(J»K) ) • F 62 

IF CJ.EO.l) GO TO A5 F 63 

CHI«J,K) - C.HI(J;K)-0T40ENS(K1+£PL/(A.*CjN)<‘0Y»*C(VB(J»K+1)+ F 64 

$ yB(J/K))+CSA(J)+(PBA(J-l/KI+?3A(J»K)l-lVB(J+l>K+X>+VB(J+l>K)>* F 65 

$ CSAt J + H*(PBA( J»K)+PSA{ J+i#K) ) ) F 96 

GO TO 50 . F 67 

C F 6R 

45 CH1(4»K) • CHIU»K)-DT*0£NS{K)*EPL/(4.*CS(4) + DY)*({VB(J>K+1)+ F 69 

$ V3(J#K))*CSA(J)*(2.*PBA(J»rO J-(VB{J+i>K+l) + V6{J+l»K))*CSA(J+l) F 70 

$ *(PBAt J»K)+P8A( J+I#K) )) F 71 

50 CO'JTIH-JE F 72 

IF tK.cl,!) GO TO 65 F 73 

IF (K.cO.KNt) GO TO 65 F 74 

00 55 F 75 

6M£5*<J) "> -DT/ (4,*PZ ) *£PL*DcNS(K> *{ ( 48 t J» K.+1)+WB( J»< ) )* (PBAtJ# F 76 



i K+l»*EPL+PBA(J»K)*EMn-(WB(J»K)+He(J,K-in*(F>aA(J»K)*EPL+PBA(J F 77 

t >K-1)*EHD) F 78 

55 CONTINUE F 79 

DO 60 4-1/JKL F 60 

60 CHI(J»K) « CHI(d#K)+GMEP( JJ F 61 

65 CONTINUE F ez 

C THICKNESS SMOOTHING F 83 

OD 70 F 84 

CHKJ^KNL) >0. F 85 

DO 70 K-1,KNL F 86 

70 PBA(J»K» » PBOUf K+1)»EPL-PB0( J^K)*£KI F 67 

DO 80 K»1>KNL F 68 

DO 75 J»3>H F 89 

75 CH1{J>K» » CHK J>K)-DT*DKYJ=(PBA{ J-2»K)-4,*PBA( J-1>K)+6.*PBA(J»K) F 90 

$ -4**PBA{J+1»K)+PBA<J+2,K))/CS< J) F 91 

CH1(2»K) • CH I(2»K)-DT*DKY<'(-3,*PBACl»K)+6«=t‘PBA(ZiK)-4,*P3A(3#K) F 92 

t +P&At4^K) )/CS (2) F 93 

CHI(1»K) « CH1(1>K)-DT+DKY*C2.*PBAJ1»K)“3.*PBA(Z>KJ+PBA<3»K))/CS F 94 

$ (1) F 95 

CHItJML»K) » CHl{JhUK)-DT+0KY*(PBA(JML-2>K)-3.*PBA(JHL-l>KJ+2*’(‘ F 96 

i PBA(JKL,KJ) /CSUHL) F 97 

60 CONTINUE F 98 

C F 99 

C UB AND VB TENDENCIES F 100 

C F 101 

DR(JHJ -0. . F 102 

DR (1) « 0. F 113 

DO iOO K-2»KNL F 1C4 

DO 95 J»2»JKL F 105 

AS - T1*UB( J»K)+T2*UB0( J#K) F 106 

AS • AS+DT+FH(J>K)+DENS(K)*C-DT/(C5A( J)**2*0Y*4, )♦( (UB(J-1/K>* F 1C7 

S CSA<J-l}+UB(J»K)*CSA(J3) + (VB(j-l,K)+CSA(J-l)+VB{4>K)+CSA(jn- F 108 

$ {UB<4#K)*CSA(J}+UeCJ+l^K)*CSAlJ+l))*(V8(J»K)*CSA(4)+VB(J+l#K)* F 109 

5 CSA(J + 1) n-DT/(4.+0Z*CSA« J> )*(<U6CJ^K)*EHI+UB{J»K-s-l)=^EPL)*CW3 F 110 

$ (J-1,K)*CS(J“1)+HB(J,K)*CSU))^(UBU>K>«£PL+UBCJ»K-1)*EHI)*(WB F Hi 

$ CJ-1»K-1)*CS(J-1)+WB( J»K-1)«CS( J3 ) J )-RAYF{K)*UB0( 4>K) F 112 

BS - Tl^VBC J^KJ+TZO'VaOC JjK) F 113 

as *> BS+DT*DENS(K )*Ue(J^K)*(UB( J-l^ K)*TNU)+UB(J + 1#K)*TNU+1) ) F 114 
$ -RAYF(K)*V30{ J^K) F 115 

C F 116 

C VB SMOOTHING F 117 

C F 115 

IF {J.EQ.2) GO TO 85 F 119 

IF (J*EQ.JHU GO TO 90 F 120 

8S « BS-DT*Di<;Y*(VB0( J-2/K)-4.*VB0U“l»K)+6.*VB0{ J#K)-4.*VB0(J+ F 121 

$ 1*K)+VB0{J+2,K))/CSAIJ> F 122 

GO TO 95 F 123 

C F 124 

85 BS • BS-DT*DKY*(+3.*VB0t2>K)-3,*VBQ{3>K)+VB0(4»K) )/CSA(23 F 125 

GO TO 95 F 126 

C F 127 

90 3S « B5-nT*DKY+{VBO(JML-2>K)-3,tVBa(JHL-l»K3+3,1'VBO(JML»Kn/ F 128 

$ CSA(JML3 F 129 

95 DR<J3 - 6S“CORlaL(J3*DT*AS F 130 

C F 131 

C SOLVE ELLIPTIC SYSTEM FOR PB F 132 

C F 133 

DO 100 J-1/4ML F 134 

RR{J»K-13 - <CHI{J»K>*E’1I*GBV(M3“CHI<O^K~ll*ePL*6flV(K-133*CNA F 135 

4 U)*(CUB(4)*DR«J)-CUB(J+1)*DR(J+1)3 F 136 

IF (K,EQ.2) RR(JjK-l) - RRtJ^ K-13-PRS (7» J3*63V(13 F 137 

100 CONTINUE F 136 

105 CALL 6LKTRI (IF LG» 1» N»AN» BN^ CN#1» JKL> AM> BM# CM> JMLj P.R> lER# WRO 3 F 139 

IFLG « IFL6+1 F 140 

IF (IFLG-13 105»105»110 F 141 

C ‘ F 142 

C COMPUTE MERIDIONAL STREAHFUNCTIQN F 143 

C F 144 

110 DO 115 K-l/KN F 145 

XBA(1>K) - 0. F 146 

XBA(JM,K) » 0. F 147 

115 CONTINUE F 143 

DO 125 J-l^M F 149 

DO x20 K»2iN F ISO 

120 XBAU+HK) . xaA(J,K)-OY+CS(J)*<CHH J»K)-(RR(J»K3*EPL-RRU»K-13* F 151 

5 EHI3)/(BVF*DT*0Z)*GBV(KJ F 152 



XBA(J-H»KyU ■ XBAt J>KNL>-DY*CS{J)*(CHI(J#KNL))/{BVF+DT*DZ)*GBV . F 153 

i tKNL) F 15A 

XBA(4+1,1) » XBA(J>H-DY*CS{J}*{CHI(J,X)-(RR{J»1)*EPL-PRS(7» J)+ F 155 

S £Hin/(BVF*DT*DI)*GEV(l) F 156 

XBA(4+1/KN» « 0. F 157 

125 CONTINUE F 158 

C F 159 

C COMPUTE NEW UB» VB» WB/ PB F 160 

C F 161 

DO 130 K-2»KNL F 162 

DO 1-30 J^SWML F 163 

AS « Tl’^^UBC J»K>+T2*UB0(J»K) F 16'« 

AS - AS+DT*F*1CJ»K)+0EK5(K)*(-DT/tC5A{J)*»2*DY+A, J*(CUB(J-1,K>* F 165 

% CSA<J-1)+UB(J»K)+CSA(J) )*(V3(J-l^K)*C3A(J-l)+VB(J>K»*CSAtJ)l- F 166 

$ (U3{J»K)*CSA(J)+UB(J+1,K)+CSA(J+1))«^(V8(J/K)*CSA(J)*VB{J+1,K>* F 167 

t CSAt J+n ) )-DT/( A.+DZ+CSAC J > )♦ t (UB( Jj ta + £rtI + UB( J#K+l)tEPU*<WB F 168 

$ (J-1»K>*CS{ J-l)+WB(J>K)*CS(J))~(UB(JfK)*£PL+UB(JjK-l)*E«n*{WB F 169 

$ (J-ljK-l)+CS{J-l)+VB(J/K-l)+CStJJ)))-RAYF(K>*U30(JfK) F 170 

UBO(J»K} « -U80( J/K)*T<t-UB(J#K)*T3 + lAS+C0R10L( J)*OT+tXBA( J»K-1 F 171 

S )*EPL-XBA(J^K)*EMI)/{DZ + CSAU) n*T5 F 172 

130 VBOCJ^K) - -VBO( JtK)*TA-T 3*VB( JjK3+T5+{XBA( J,K-1)+EPL-X8A{J»K)+EHI F 173 

S)/(OZ«CSA(Jn F 174 

DO 135 J»1jJ1U F 175 

DO 135 K“1»KNL F 176 

135 WBO(JjK) “ -WB0CJ#K)*T4-T3*WB(4^K)+T5*(XBA( J#K3-XBA(J+1»K))/(0Y«CS F 177 

$(J)) F 176 

DO 140 J«1^JHL F 179 

PBA(J»i) - PRS(7^J) F leo 

UB0(J»1) * PRSt5»JJ F 161 

00 140 K«2,KNL F 162 

140 PBA(J>K) » RRTJ»K-1) F 183 

DO 145 K«lfKN F 164 

145 PBAtJHfK) • PBA{JKL/K> F 165 

DO 150 J=l/4tt F 186 

PBAtJ^KN) • PBAU,KNL)*EMI/EPL F 167 

VBO(JfKN) - VBO( J/KNU+EKI/EPL F 165 

150 UB0(4,KJ<) • U30( J#KNL)*EHI/£PL F 169 

DO 155 F 190 


DO 155 K"1»KN 

CHI(J»K) - T51'PBA(J»K)-T4*PBQ( J^K)-T34PB( J>K) 
PBO(J#K) - PB(J,K) 

P8( J»K) •= CHI(J#K) 

CHKJjKJ - WBOCJ^K) 

WB0(4»K.) » WBtJjKJ 
WB(4/K) - CHItJ/K) 

IF (J.EQ.l) GO TO 155 
CHKJ^K) • UBO(J>K) 

U30(J»K) « UB(4>K) 

U8(J>K) •'-CHK4/K) 

CHICJ^K) • VBOtJ^K) 

VBO(J>K.) • VS(J»K) 

VB(J»KJ « CHI(J>Kl 
155 UB(J^K) » UB(J»K)*DENS(K) 

RETURN 

END 


F 1«1 
F 192 
F 193 
F 194 
F 195 
F 196 
F 197 
F 198 
F 199 
F 200 
F 2C1 
F 202 
F 203 
F 2C4 
F 2o5 
F 206 
F 2o7- 



SUBROUTIKE EDDY BB*CC> ALB, 0 VF# ANGj CI» CS> CS A» CaR»DC A# DEL, 0£NS> 
$DKY, DTOOY,DT, DZ, DY, EHI, EPL,GAH, GH, ICT, IFLG, IFD, JH,KL'L,K,S,t^GC, NU, 
SRAYF, PLl, P, Q,S’J,SV,TU#TY,TNA,TN8#XH,X«A, PB,PHIA»PHl,PHIO,R, UB,U,UO 
$,V>VO»H,WRO,CDUM» AN» BN, CN,G8 V, S Yh, Q5, KZ ) 

DIHENSIOH AA{1), B&d), CC(1J, WRO(l), PLIUI, PHIA{JH,1), 

S PHIO{JH,X>, PHI(JM,1), U<JK,1), V{JH,1>, VD(JM,1), ANG 

$(1J, XMA(l), UB<J«,1), CS(1>, CSA(l), TNA(l), PCI), Q(l>, CORll), 
SDELUI, MGCdJ, NL)tl>, DCA(XJ, TU(1), SU(L), TVU 

$), SV(1), ALBU)# OENSd), 6AH(1), RAYFdl, TNBCi)/ X/ICl) 

$, CDUH{JH,1>, AN(1>, BN{1), CNdl, GBVd), QSUH,1) 

REAL NGC,KZ 
INTEGER S,SYH 

COMPLEX BB,R, PH I A, PHI 0, PHI, U, UD, V, V0,CI,P,0, NU, Gft, TU,TV, AV, BV, PLl, 
$W#CDUH,QS 
KN • KNL+1 
N e KNL-1 
JHL • H+i 
DO ID 

DO 10 K-2,KN 

10 UB(J»K) » UB(J,K1 /DENS(K) 

C CHOOSE LEAPFROG OR FORWARD DIFFERENCE 

ICT » ICT+i 

IF (ICT.LT.IFD) GO TD 15 
T1 - T3 - 1. 

T2 * TA B 0* 

T5 « 2. 

GO TO 20 
C 

15 Ti • T2 • .5 
T3 » 2. 

TA - 1. 

T5 - 

20 CONTINUE 
ISW « 1 

C THICKNESS TENDENCY 

DO 30 J-1,JH 
DO 25 K-1,KNL 

25 PHIA(J,K1 » PHIQ( J»K+l)+tPL-PHIO(J,K)*EHI 

30 PHIA(J,KN) » «0.,0.) 

DO 35 K^lrKNL 
00 35 J-2,Jf'L 

CDU>1{J,KJ “ (T2“0CA(K))*PH1A( J,KJ + (Tl-OcNStK> + EPL*Grt{Jl + JUBCJ, 
$ K+1) +U8( J, K) ) )^(PHH J,K+l)*cPL-PHIL J, K) ♦£«! )-DTODY*OENS (K)*EPL 

i *(V( J-1,K+1)+V( J,K+1 )+V(J-l,K*+V( J,K) )*C(P3( J-1,K+1)-PB(J,K+1) 

S )f'£PL-(P8 (J-1,K)-PBU,K)> + EHIJ+DT*Q5( J,K) 

IF (K.EQ.l) GO TO 35 

CDUHU,K) - CDUH( J,Kl“DT*D£NSLK>*EPL/(A.»DZ)*W( J,K)*(((PB( J~l, 
$ K+2)t-PB( J,K+2) )‘tEPL-(PB(J-l,K + l) + P3 (J,K+i) > *EHI ) * = PL**2- ( ( P B ( J 

$ “1,K) + PBJ J,K))*EPL-(PB{ J-l,K-U + PB(J,K-i) J*ErtI)*EMI**2) 

35 CONTINUE 

C THICKNESS SMOOTHING 

DO 45 J»2,JK 
00 40 K-3,N 

40 CDUH(J,K) ■ CDUH(J,K)-DT*KZ+(PHIA( J,K-2 J-4.+PHIA(J,K-l)+6.*PHIA 

t ( J,KJ-4,*PHIA( J,K+1)+PHIA(J,K+2J ) 

CDUK( J,2) « COUM( J,2»-DT*KZ*(-PHIA( J#l)*3**PHIA(J,2)-3,*PrtIA'{J,3 
$ )*PHIACJ,4)> 

COUH(J,KNL> - CDUHt J,KNLJ-DT+KZ*(PKIA( J,N-1J-3.*PHIA(J,N)+PHIA( J 
$ ,KNL)*3.) 

45 CONTINUE 

DU 55 K»1,KNL 
DO 50 J«3,M 

50 COUH(J,K) « CDJM( J, K)-DT’tDKY<'(PHlA(J-2,K)-4.*PHIMJ-i,K)+6.*PHIA 

S (J,K)-4.*PHIA( J+1,K)+PHIA( J+2,K) )/CSA{ J) 

CDUH{2,K) • CDUHL2, K) -OT*DKY* ( 3,»PHIA (2, K)-3 .+PH1 A{ 3,K )+PriI A (4, K 
S n/CSA(2) 

CD'JH(JHL,K) « CDUPt Ji1L,K )-0T*0KY*{3.*PHlA( JNL,K)-3.*PHIA<JHL-1,K 
J )+PHIA(JML-2,K))/CSA(JNL» 

55 CONTINUE 
C MOHcNTUH TENDENCY 

60 DO 115 K«1,KNL 

IF ( (K.EQ.l).AMO. (ISW.EO.l) J GO TO 115 
DO 110 

ANGCJJ ■ (XHA( J1*U3(J,K)+XMA(J+1)»U3{ Jtl,K) )/2. 

AV « T1*U(J,K»+(T2-RAYF fK) )*UDC J,KJ-DT*0£NS(K)* (CI*ANG{J)*U(J, 
$ KH-V( J,K)/(CS( J)*DY) + (U3 (J,N )*CSA ( J J-U6 ( *CSA ( J+1 ) )+ ( W < J + 


G 1 
G 2 
G 3 
G. 4 
G 5 
6 6 
G 7 
G 8 
6 9 

G 10 
G - 11 
6 12 
G 13 
G 14 
G 15 
G 16 
G 17 
6 18 
G 19 
G 20 
G 21 
G 22 
G 23 
G 24 
G 25 
G 26 
G 27 
G 28 
6 29 

G 30 
G 31 
G 32 
G 33 
G 34 
6 35 

G 36 
G 37 
G 3P 
G 39 
G 40 
G 41 
G 42 
G 43 
G 44 
G 45 
G 46 
G 47 
G 46 
G 49 
6 50 

G 51 
G 52 
G 53 
G 54 
G 55 
G 56 
G 57 
G £8 
G 59 
G 60 
G 6l 
G 62 
6 63 

G 64 
6 65 

6 66 
6 67 

G 68 
6 69 

G 70 
G 71 
G 72 
G 73 
G 74 
G 75 
6 76 



» 1»K)*(UB(J+1*K-H)*£PL-UB{J+1>K)*EHIJ+M<J/KJ*{UBU#K+1I*EPL--U8 6 77 

$ <J»K)+EHI)>/(A,*D2)) 6 78 

IF (K.GT.l) AV - AV-DT*0cNS(K)+(W(4+l#K~i)*(UB(J+l>K)*EPL~UBU 6 79 

$ +1/K-1J*EHI)+W( J>K-l)*^(U8lJ/K)«cPL»U3{J#K-l)*EKI) )/(4.+ DE) G 60 

BV “ Tl*V{J/K)+tT2-RAYFCK)}*V0(J/KJ-DT*0ENS(K>*{CI*ANG< J)*V(d# G 81 

$ K)-U(J#K)*<UB(J/K)*TNA{J)+UfiCJ-H,K)*TNB{ Jin 6 62 

C HOHENTOft SHOGTHING G 63 

IF (K.EQ.il GD TO 75 G 64 

IF (K.EQ,2) GO TO 65 G 65 

IF (K.EO.KNL) 60 TO 70 G 86 

AV « AV-DT*KZ*<UD(J»K-2)-4.*UD( J>K-1)+6.*UO(J»K)-4»*UO(J/K+1)+ 6 87 

$ Ua{J,K+2)) G £8 

BV - BV-DT«^KZ*( VO(Jj.K-2)-4**VO{J/K-l)+6,*VO( J»K)-4,*V0(J>K+1)+ G 69 

$ V0{J,K+2)) G 90 

60 TO 75 G 91 

C 6 92 

65 AV » AV-DT*KZ«(“U0(J,1)+3.*UQ{J,2I-3.*U0< Jj3)+UQ{ J^4U G 93 

BV » BV-DT*KZ*(-V0(J/11+3.*VQ(J>2)“3.*V0( J/3)+V0{ J/9) ) G 94 

GO TO 75 6 95 

C G 96 

70 AV - AV-DT*KZ*(U0(J#N-1)-3.*U0( J^N)+3.*U0{ J»KNU) G 97 

BV ' BV-DT*KZ*(V0(J»S-1)-3.*V0( J»NJ+3,*V0(4»KNL)) 6 98 

75 IF UoEQ.13 GO TO 80 6 99 

IF (J.EQ.2) GO TO 85 G 100 

IF (J.EQ.JML-1) GO TO 90 6 101 

IF U.EQ.JHLl GO TO 95 G 102 

AV - AV-DT*DKY*(U{HJ-2,K)-4.*UO(J-l,K>+6.+UOU>K>-4.*UOU+l»K) G 103 

$ -+U0( J+2>K))/CS( J) G 104 

BV - BV-DT*DKY*( V0( J-2#K»-4.*VO(J-1»K)+6.*VOIJ/K)-4.*VO( J+l/K) G 105 

$ +V0( J+2#K) )/CS(J > G 106 

GO TO 100 6 107 

C 6 108 

80 IF (S.EO.l) CF * 2.S IF (S*NE.l) CF « 4. G 109 

AV - AV>DT*DKY*(CF*U0(l»K)-3.*U0(2#K)+U0(3^Kn/CSa) G 110 

BV - 6V--DT*DKY*(CF*VO£l,K)-3, + VD(2>K)+VO(3,Kn/CSa) 6 111 

GD TO lOo 6 112 

C G 113 

65 IF (S.EQ.IJ CF » 3.S IF (S.NE.l) CF » 5. G 114. 

AV •" AV-DT*Di4Y*(-CFfUO{i,K)+6,*=UQ(2»K>-44*UO(3»K)+UO(4»Kn/CS G 115 

J (2) 6 116 

BV » BV*-DT»DKY+(-CF + VO(l#K)+6.*VO{2»K>-4.*VO(3/K)+V£JC4»K)) /CS G 117 

i (2) 6 118 

GO TO 100 G 119 

C G 120 

90 IF (S.EO.l) CF - 3.S IF (S.NE.l) CF - 5. G 121 

IF (5Y«.E0.<.) CF • 5, G 122 

AV « AV-0T*DKY*(-CF*UD(JHLjK)+6.+U0( JKL~1^K)-4.*UQ(JHL-2»K)+U0 G 123 

$ (JHL-3#K))/CS(JHL-1) G 124 

IF (SYM.FQ.l) CF - 3. G 125 

6V - BV-DT*D<Y*(“CF*VOUKL>K)+6.*VOUKL-1»K)-4,*VOUHL-2»K)+VO G 126 

S ( JHL-3>K) )/CS(JHL~l » G J.27 

60 TO 100 G 12S 

C G 129 

95 IF (S.EO.l) CF » 2.S IF (S.NE.l) CF - 4, G 130 

IF (SYH.EQ.n CF - 4, G 131 

AV - AV-DT*0KY*(CF*U0{JHL^K)-3.*UD( Jf1L-i#K)+U0(JML-2»K))/CS G 132 

$ (JML) G 153 

IF (SYH.EO.IJ CF - 2. G 134 

aV - BV-DT*DKY*(CF*VOI4rtL»K)-3.*VO(JHL-l#K)+VO<JML-2»K> l/CS G 135 

$ (JNU G 136 

100 CONTINUE G 137 

C COMPUTE FORCING TERM IN ELLIPTIC EQUATION 6 138 

IF USW.E0.2) GO TO 105 G 139 

?(J) - AV+COR( J)*DT*3V G 140 

0(4) * BV“CaR£4)*DT*AV 6 141 

IF (J.EO.l) GO TO 110 G 192 

R(J-1#K“1> - COUM< J#K)*EMI*GBVtK)“CDUfUJ>K.“l)*EPL*GaV(K-l>+DcL G 193 

$ {J)*{(NGC(J-1)*0( J-'1)-NGC(J)*Q(J) )/OY+NU(J-l)*P(J-l)+NU(J)*P(J G 144 

4 )) G 145 

60 TO 110 G 146 

C KEW VALUES FOR U AND V G 147 

105 UO(J»K) - -UOC J^K »*T4+T5*(-TU( J)*(PHIA(J+1>K)+PHIA( J>K))-SU(J) G148 

4 ■ *{PHIA(J.N-)-.»riIA(4*l»'<ntAV*‘GAfl{JJ+ALb{J)»5V)-T3*U( J/K) G 149 

VOCJ^O . -VD(j.K)*T 44T5*^(TV (4 ) + (PNlA(J + l,K)+PHIA(J>K))-SV(J)* G 150 

$ (PHIA(J»K)-PHIA(4+l>K))+BW*GAh(J)-ALB{J)*AV)-T3*V(J>K) G 151 

110 CONTINUE G 



115 


120 

C 

125 

130 

135 

140 

C 

145 

150 

C 


155 

160 


CaHTIKUE 

IF aSW.EQ*2) GO TO 145 
DO 120 J»1>K 

PHIA(J+1/KNL+1) « (0.»0.) 

R(J»11 - R( J#i)-PLl(0*l>«G6Vm 
INVERT ELLIPTIC EQUATION FOR PHI 

CALL CBLKTRI CIFLG^ 1> AN^ BN# CN« 1> H^AA^ CC> IER> WRO ) 

IFLG • IFLG-5-1 
IF CIFLG-1) 125#125#130 
CONTINUE 
DO 135 J»2>J«L 

PHIA(J»ll » PLl(J) 

DO 135 K*2#KNL 
PHIA(J»K) « ft(J-l#K-l) 

DO 140 K-1#KNL 

PHIA(i/K) * (0i#0.) 

PHIA(JN>KJ » C0.#0.) 

ISM - 2 
GO TO 60 

CONTINUE 
DO 150 J-2#JHL 
DO 150 K-1#KNL 

H( J#K) » (C0UH(J#K)-(PHIA(J#K+1)=*EPL-PH1A(J#K)*EHI) )/{ 8VF+DT*D21* 
$G8V(K) 

NEW VALUES FOR DEPENDENT VARIABLES 
DO 155 

DO 155 K-1>KN 

PHIA(J#K) » T5*-PHIA( J>KJ-T4*PH10<J#KJ-T3*PHICJ#K) 

PHIO(J#K) • PHI(J#K) 

PHI(J>K) » PHIA(J»K) 

PHIA(J#K) - UDU»K> 

UO(J»K) « U(J#K) 

U(J»K) »« PHIA(J»K) 

PHIA<J#K) » VQtJ#K) 

VO(J#K) - V(J>K) 

V(J#K) • PHIA{J#K) 

CONTINUE 

DO 160 K«1#K.N 

- -U(J«L»K) 

V{JM#K) » V(JHL#K) 

00 160 J«>2»JH 

UBtJ/K) “ Ua(J# K>*D£NS(K) 

CONTINUE 

RETURN 

END 
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ORIGSNAURAGE ■ IS ' 



SUBROUTINE FLUX (JK> KN/DY» D2>CSA»CS» DENS# £PL» E«I, U/ Vj W> P» T#F M, FT, KZ H 1 

i,UBO, VT,CNB) H Z 

DIKENSION CSA(l), CS(1), DENSd), U(JH,1), W(JN,1), P(JM> H 3 

$1), T(JH,1), FT(JM,X), H A 

REAL KZ H 5 

DIMENSION VT(JK,1), CNB(l) H 6 

COMPLEX U>V,VI,P,T H 7 

KNL « KN-l H 8 

OML - JH-1 H 9 

N ■ KNL-1 H 10 

DO 10 H 11 

DO 10 K«1,KNL H 12 

10 T(J,K) » P(d,K+l)*EPL-P CJ,K)*EMI H 13 

DO 15 K»1,KNL H lA 

VT(1,K) - VT(JH,K) * 0. H 15 

DO 15 H 16 

15 VT(J,K) * 0£NS(K)*£PL«^CSAfJ) + (R£ALCV( J-1,K+1)4V(J,K^1^+V{J-1,K)+V H 17 

S(J,K) KREAL (TCJ>K) 1 +A1KAG ( V ( J-1, K41 )+ V ( J, K+i J + V( J-1, K ) +V ( J, K ))* H 18 

SA1KAG{T{J,K) ) >/2. H 19 

DO 20 J-1,JML H 20 

DO 20 K»1,H H 21 

20 FT(J,K> - ••CN8{ J )t=(VTtJ,K)-VT( J+l,Kn H 22 

DO 25 J»1,JHL H 23 

VT(J,KN) “0. H 2A 

VT(J»KNL) « 0. H 25 

DO 25 K-1,N H 26 

.25 VT(J,K) - REAL{H<J4l,K>>*PEAL(T(J*l,K))+ft£AL(W<J,K))*REAL(T(J,Kn+ H 27 

$AIHAG(rf(J+l,K))4AIHAGtT{J+l,K) }+AIMAG(W( J,KJ ) +AIHA6 (T ( J, KU H 2B 

DO 30 K»2,N H 29 

DO 30 J»1,JHL H 30 

30 FT(J,K) « FT(d,K)-DENS(K»*EPL’P(VT{J,K+l)-VT( J,K-in/(2. + DZ> H 31 

DO 45 J»2,JM H 32 

DO 45 K-l,KNl H 33 

FM(J,K) « -DENSIK)*UREAL(U( J-l,Kn+RtAL(Vtd-l,K) } + AIMAG(U(J-l H 34 

$ ,K))*AIMAGl V(J-1,K)))*CS(J-1)**2-(R£AL(U(J,K))4REAL(VU,K)}+ H 35 

S AIHAG(U(J,Kn*AIHAG{V«J,K)))*CS(JJ**2)/<CSA(J)**2*0Y>+2. H 36 

IF (K.EQ.l) GO TO 45 H 37 

FHU,K) - FM( J,K)-DENS{K)*t tREAUU(J-l,i<)+U{J,Kl )*EMI+REAL{U< J H 38 

$ -1,K+1)+U<J,K+1) ) *EPL )4REAL(W< J,K))+( AIMAG(U(J-1,K)4U{J,K) )* H 3° 

S £MI + AIHAG<U( J-1,K+1 )+U{ J,K+1) J»CPUfAlMAG(HU,K) )-(REAL(U( J-1, H 40 

$ K)+U(J,K) )*EPL+RFAL (U(J>-l,K-l)+U(J,K-in4EMn*R£AL(rftJ,K-in- H 41 

S UIMAGa)(J-l,K)+U(J,Kn*EPL + AlHAG(UlJ-i,K-l)+U(J»K-l))*EHI)+ K 42 

S AIMAGCVK n/(2.*DZ) H 43 

IF <K.EQ,2> GO TO 35 H 44 

IF (K.EQ.KNLI GO TD 40 H 45 

FM(J,K) - FM (J,K)-KZ*(US0(J,K-2)-4.*UB0(J,K-l)+6t'»UB0<J,R»-4,+ H 46 

$ UBO( J,K+1)+UB0< J, K+2n H 47 

GO TO 45 H 48 

C H 49 

35 FM(J,2) « FM(J»2J-KZ*(“U30{J,l)+3,*UB0<J,2)-3.*UBO{J,3)+UBD{J, H £0 

i 4)) H 51 

GO TO 45 H 52 

C H 53 

40 FH{J,KNU « FHU,KNL)-KZ*(UBDU,KNL-2)~3.*UB0(J,KNL-1)+3.*U30 H 54 

$ (J,KML)-UBO( J,KN) ) H 55 

45 CONTINUE H 56 

RETURN H 57 

END H 58- 



SUBROUTIME AOUT(TIHE^DT>KN# JM>DaH»D£NS»IPHAS»i#JHL>DZ>RHZ»CaR»DY# 
SCI,XH1 /IENOjITIKE»«T, IMAT^P l^ PB»U1,UB»V1# VB^W1>WB»FH,FT/S»0S#CS>QB 
$>CZT»SZT,COT»TQ) 

DIHcNSION OUMUtili OENS(KN)# H9)f PB(JMjKN), U8(JM>KN)> VB(JM^KN) 
%f WB(JH,KN)# CS(JKJ» QB<JH#KN># C2T(JK#Kt<J, SZT(JH>KM)# CDT(JH>KN) 
S> TQ(JML#16), TDQ{18/16), LAT(19> 

DIHEMSION IPHAS(JH>» COR(J«>, XHKJMl# F.KJM^KN)# FKJH^KN) 

COMPLEX Cl, P1(JH.KN),U1(JH»KN)»V1{JH>KN)# WH JM# 16 >, QS < JH, 16 ) 

INTEGER S 

PI • 2.*ASIN(1.) 

DAY - (TIME )/(3600.*2^«) 

KNL e KN-1 
WRITE (6,105) 

WRITE (6,130) DAY 
DO 15 KK“1,KN 
K • KN-KK+1 
DO 10 J»1,JH 
DU.M(J) » 0,0 
DUM(J) • VB{J,K)*DENS(K) 

10 CDRTIHUE 

22 * Z(K)+16000, 

WRITE (6,190) 

WRITE (6,110) ZZr (DUH( J), J-1,JM) 

15 CONTINUE 

WRITE (6,105) 

WRITE (6,135) 

DO 30 KK-1,KNL 
K « KN-KK 
DO 20 

20 DUM(J) « RH2*(PB ( J,K+1)*DENS<K+1>-PB{J,K)*0ENS(K) )/0Z 

ZT « Z(K)+DZ/2, +16000, 

SUN - 0# 

DO 25 J-l,JMl 

TDQ(J,K) » DUH(J)-TQ( J,K) 

25 SUN « 5UN+DUH(J) *CS(J) 

SUN « SUH+PIY36. 

WRITE (6,190) 

WRITE (6,140) ZT, (OUK(J),J-lr JML),SUN 
30 CONTINUE 

WRITE (6,105) 

WRITE (6,145) 

DO 35 KK=1,KN 
K • KH-KK+1 
ZZ « Z(K)+16000, 

WRITE (6,190) 

WRITS (6,125) ZZ, {U8 (-J,K)> J-1, JH > 

35 CONTINUE 

DO 40 L-1,19 
LL - (lO-D+lO 
LAT(L) • LL 

40 CONTINUE 

WRITE (6,150) LAT 
WRITE (6,105} 

ViPITE (6,155) 

DO 50 KK»l,KN 
K - KH-KK+1 
DO 45 

45 OUM(J> - WB(J,K)*DENS(K)41.E3 

ZZ • Z(K)+16000.+DZ/2, 

WRITE (6,190) 

WRITE (6,120) ZZ, (DUH(J ), J«1,JML> 

50 CONTINUE 

WRITE (6,105) 

WRITE (6,160) 

DO 60 KX»1,KN 
K « KN-KK+1 
DO 55 J»1,JH 

55 DUM(J) • CB(J,K>*0ENSCK)*86400.*RHZ/0Z 

ZO • Z( K ) +DZ/ 2,+lo vwO, 

WRITE (6,190) 

WRITE (6,115) ZD, (DU«(J ), J»1,JHLI 
60 CONTINUE 

WRITE (6,10f) 

WRITE (6,100) 

DO 65 KK«1,KN 
K ■ KN-KK+1 


I 1 
I 2 
I 3 
1 4 

I 5 
1 6 
I 7 
I 8 
1 9 

I 10 
I 11 
I IZ 
I 13 
I 14 
I 15 
I 16 
I 17 
I IB 
I 19 
I 20 
1 21 
I 22 
I 23 
I 24 
I 25 
I 26 
I 27 
I 28 
I 29 
I 30 
I 31 
I 32 
I 33 
1 34 

I 35 
I 36 
I 37 
1 38 

I 29 
I 40 
I 41 
I 42 
1 43 

I 44 
I 45 
I 46 
I 47 
I 48 
I 49 
I fO 
1 51 

1 52 

I 53 
I 54 
I 55 
I 56 
I 57 
I 58 
I 59 
I to 
I 61 
1 62 
I 63 
I 64 
I 65 
I 66 
I 67 
I 66 
I 69 
I 70 
I 71 
I 72 
I 73 
I 74 
I 75 
I 76 



ZO ■ 2(K>+DZ/2. +16000, I 77 

WRITE (6>190) 1 76 

WRITE (6/115) ZD, (SZT(J/K)»J«*1/ JKL) I 79 

65 CONTINUE I 60 

WRITE (6,lo5) I 61 

WRITE (6,165) DAY,S 1 62 

DO 75 KK-1,KN I 63 

K • KN-KK+1 I 84 

DO 70 I 85 

OUH(J) » CABS(P1( J,K) ) I 66 

IF (DUMIJI.tO.O.) GO TO 70 I 87. 

DUH(J) « DUH( J)+DEMJ{K) /9,6+2, I 68 

IPHAS(J) « ATAN2( AI«AG(P1(J,K) ),REAL(P1(J,K) ) )*160./PI 1 89 

70 CONTINUE I 90 

ZD « Z(K.)+16000, • 1 91 

WRITS (6,160) ZO, (DUM( J ),J-1,JHL) I 92 

WRITE (6,185) UPHAS(J),J«1,JML) I 93 

WRITE (6,190) I 94 

75 CONTINUE 1 95 

WRITE (6,105) I 96 

WRITE (6,170) I 97 

DO 85 KK«L,KNL I 96 

K • KN-KK I 99 

DO eO I 100 

60 DUH(J) FM( J,K)«^0ENS(K)+1.E6 I 101 

DUM(19) » 0. I 102 

ZD » Z(K)+16000; I 103 

WRITE (6,190) I 104 

85 WRITE (6,125) ZD, (DUN(J), J»l, JH) I 105 

WRITE (6,105) I 106 

WRITE (6,175) I 107 

DO 95 KK-1,KNL 1 106 

K « KN-KK I 109 

00 90 I no 

90 OUM(J) - FT(J,K)*DLNS (K )<=RKZ/DZ+1.E6 I 111 

ZD “ Z(K)+16000. I 112 

WRITE (6,190) I 113 

95 WRITE (6,125) ZO, ( DUH( J ), J»l, Jf.) I 114 

KT « 1 I 115 

IMAT “0 I 11b 

RETURN I 117 

C I 116 

C . I 119 

100 FORMAT (+ SOLAR HEATING 6Y OZONE f) I 120 

105 FORMAT (IHl) I l^i 

110 FORMAT <-3PF6.1, 19(0PF6.2) ) I 122 

115 FORMAT (-3PF6»1,18(0PF6»2) ) I 123 

120 FORMAT (-3PF6.1 , 16 (0PF6. 1 ) ) I 124 

125 FORMAT (-3PF6 ,1, 19(0PF6 ,i ) ) I 125 

130 FORMAT (29H MEAN MEPlDiaNAL WIND DAY« ,F6,2) I 126 

135 FORMAT (20H MEAN TEMPERATURE ) I 127 

140 FORMAT ( -3 PF6 . i, 16 ( OP F6 . 1 ) , 6X , FIO. 4 ) I 128 

145 FORMAT (IftH KEAN ZONAL WIND ) I 129 

150 FORMAT (/,6X,19I6) I 130 

155 FORMAT (26H VERTICAL VELOCITY, KH/S ) I 131 

160 FORMAT (30H RADIATIVE HEATING K/D ) I 132 

165 FORMAT (Z3H GEQPDTHNTIAL, DAY • ,F6.2,12H WAVENUMBER ,14) I 133 

170 FORMAT (30H EDDY MOMENTUM FLUX DIVERGENCE) I 134 

175 FORMAT (2faH EDDY HEAT FLUX DIVtRGENCE) I 135 

180 FORMAT (-3PF6,!, 19 (0PF6.1 ) ) 1 136 

185 FORMAT (6X,19I6) I 137 

190 FORMAT (IH ) I 156 

end I 139- 



